Equations of State with Group Contribution Binary Interaction Parameters for Calculation of Two-Phase Envelopes for Synthetic and Real Natural Gas Mixtures with Heavy Fractions

— Three equations of state with a group contribution model for binary interaction parameters were employed to calculate the vapor-liquid equilibria of synthetic and real natural gas mixtures with heavy fractions. In order to estimate the binary interaction parameters, critical temperatures, critical pressures and acentric factors of binary constituents of the mixture are required. The binary interaction parameter model also accounts for temperature. To perform phase equilibrium calculations, the heavy fractions were ﬁ rst discretized into 12 Single Carbon Numbers (SCN) using generalized molecular weights. Then, using the generalized molecular weights and speci ﬁ c gravities, the SCN were characterized. Afterwards, phase equilibrium calculations were performed employing a set of ( nc þ 1) equations where nc stands for the number of known components plus 12 SCN. The equations were solved iteratively using Newton ’ s method. Predictions indicate that the use of binary interaction parameters for highly sour natural gas mixtures is quite important and must not be avoided. For sweet natural gas mixtures, the use of binary interaction parameters is less remarkable, however.


INTRODUCTION
Information provided by vapor-liquid equilibrium is important in gas processing [1,2] and the related environmental protection [3]. While phase envelope calculations for synthetic and lean natural gas mixtures with a few known components do not often pose problem [4], for real gas mixtures with unknown fractions calculations are increasingly more complicated [5]. Often the presence of a third phase cannot be excluded a priori [6] and if water is present, water dew point curve together with hydrate phase calculations need to be simultaneously addressed.
A real gas mixture collected from well-head is often reported by a few known components plus an unknown fraction [7]. Prior to phase equilibrium calculations, the unknown fraction needs to be discretized into a certain number of Single Carbon Numbers (SCN) and then the SCN need to be characterized by available methods [8][9][10]. With the advances of molecular-type equations of state, characterization of petroleum fluids has been remained an active field of research [11,12]. As the equilibrium measurements for real gas mixtures are difficult to measure, expensive and limited, numerical methods using Equations of State (EoS) have become imminent [13,14].
Use of reduction methods [15,16] and change of the space variable from temperature/pressure to density are just a few [17,18].
Numerical methods take the advantages of EoS to describe vapor and liquid phases. The field of EoS has been improving since the seminal work of van der Waals [19]. Lately, EoS with Helmholtz free energy form found particular attention in natural gas equilibrium calculations [20]. The European Gas Research Group (GERG) developed GERG-2008 as a standard tool for lean natural gas property calculations [21]. In terms of accuracy, comparisons made with other EoS [22,23] did not rule out the use of simpler EoS in phase envelope calculations. Especially, when real natural gas mixtures with multitude of constituents are encountered, the use of cubic EoS is still the most convenient. As such, attentions have been made to improve the Pressure-Volume-Temperature (PVT) of the cubic EoS [24,25]. In particular, more attentions have been made to make established EoS such as Peng-Robinson (PR) [26] or Soave-Redlich-Kwong (SRK) [27] versatile and more accurate. EoS like Predictive SRK (PSRK) [28], Volume Translated PR (VTPR) [29], Universal Mixing Rule for PR with UNIFAC (UMR-PRU) [30], and Predictive PR (PPR78) [31][32][33][34][35][36][37] are just a few examples. The latter, however, found especial attention in hydrocarbon systems as it incorporates a predictive temperature-dependent group contribution binary interaction parameter model into the PPR78 EoS [38]. The ease of use, satisfying accuracy and potential wide spread applications to a wide variety of systems containing paraffines, naphthenics, aromatics, nitrogen, sulfur compounds and carbon dioxide makes PPR78 attractive in gas industry. Recent applications of PPR78 extended the use of PPR78 with reduced flash methods [39]. Additionally, with minor adjustment the group contribution binary interaction model applied to PPR78 can be coupled with similar EoS but different PVT relationships [36,40].
Binary interaction parameters effectively improve the EoS power in describing the phase equilibria of systems containing carbon dioxide, hydrogen sulfide, asymmetric hydrocarbon mixtures, and systems containing associating fluids [41][42][43][44]. For light natural gas mixtures [45] devoid of appreciable amount of sour gases, the use of binary interaction parameters is trivial as pointed out by Nasrifar et al. [4] and quite recently by Skylogianni et al. [46]. Mørch et al. [47] did not use any binary interaction parameters and obtained satisfying results in modeling the dew points of light synthetic natural gas mixtures.
Real natural gas mixtures, however, contain multitude of components varying from paraffins to naphthenics and aromatics [7]. Natural gas mixtures may contain appreciable amounts of nitrogen, carbon dioxide and sulfur compounds as well. In past, describing such complicated mixtures were being accomplished by incorporating temperature-indepen-dent binary interaction parameters between the few light components and/or use of Chueh and Prausnitz model [48] between the heavy hydrocarbons including SCN [7]. The latter only takes the benefits of critical volume for hydrocarbons and its accuracy is often taken for granted. As such, PPR78 finds applications in hydrocarbon processing knowing that PPR78 uses a group contribution model for estimating binary interaction parameters and that it can be applied to almost all structurally known components present in real natural gas mixtures. Nevertheless, real natural gas mixtures with unknown fractions contain SCN with unknown structure.
In this work, we extend the application of PPR78 in calculating phase envelopes to real natural gas mixtures with unknown fractions. The binary interaction parameter model undergoes stringent test for temperature dependence and presence of long chain molecules. Then, with some minor adjustment, as pointed out by Jaubert and co-worker [36,40], the group contribution model will be used with the Nasrifar and Bolland (NB) EoS [49] and the SRK EoS. The coupling of SRK EoS with the binary interaction model will be called PR2SRK here on after Jaubert and Privat [36]. These two EoS are also used in calculating phase envelopes as well.
Once the EoS, mixing rules, the group contribution model and the characterization method are presented, the calculation method [6] will be provided and comparisons will be made with experimental data. The strengths of the methods are elucidated and the paper is concluded with specific remarks.

EQUATIONS OF STATE
The PVT of the EoS studied in this work can be expressed by where p is the pressure, R the gas constant, T the temperature, v the molar volume, b the molar co-volume, a the attraction parameter and d 1 and d 2 are the two specific constants that vary depending upon the EoS. The values of d 1 and d 2 for PPR78, PR2SRK and NB EoS are presented in Table 1. The values of a and b are determined by critical temperature, critical pressure and acentric factor. While b is temperature independent, a is a similar function of temperature Specific constants for the studied EoS [31,36,49].
for the three EoS. With the exception of NB whose temperature dependence was augmented for supercritical region [50], the subcritical temperature dependence of the three EoS has the same form. The expressions for a and b can be found in [31,36,49], respectively. The van der Waals mixing rules were employed to extend the EoS to the mixtures: where x is the vapor or liquid mole fraction, and k ij the Binary Interaction Parameter (BIP). Jaubert and co-workers [31][32][33][34][35] developed the following temperature-dependent, group contribution BIP model for the PPR78: where indexes l and k run over the groups, ng is the number of groups, a ik is the fraction of molecule i occupied by group k and A kl and B kl are the interaction parameters between groups k and l. The interaction parameters between groups were found in a series of work [31][32][33][34][35] and can be found in a complete form in [51]. While Equation (4) was derived for the PR78 EoS, the authors outlined how to apply it to the NB [48] and SRK EoS [27]. It is just enough to multiply the group interaction parameters with a constant À that is, where h is the multiplier of the proposed EoS. The values of h are provided in Table 2.

METHOD OF CALCULATIONS
Phase envelope calculations were described in details by Michelsen [6]. A brief opted discussion of the method follows. Phase boundaries for vapor-liquid equilibria are calculated from the mole fraction constraint: where z j is the mixture mole fraction and K j the equilibrium ratio. K j is defined as the ratio of stable phase to incipient phase (vapor or liquid) mole fraction. In terms of fugacity coefficients, K j can be expressed by: where the superscript incip stands for the incipient phase and nc for the number of components. Equations (7) and (8) form a complete set of (nc þ 1) equations. Specifying T (or p), one only need to solve the above system of (nc þ 1) equations for K j 's and p (or T), respectively. As such, we define the solution vector by: where Y represents lnT (or lnp) once specifying p (or T), respectively. The vector of function values is expressed by: A particular solution of Equation (10) is found by Newton's method, i.e., J ðmÞ Da þ g ðmÞ ¼ 0; ð11Þ where m refers to the iteration m, J is the Jacobian matrix calculated at a (m) and the Jacobian elements are J ji = ∂ g j /∂ a i . The Jacobian elements were calculated analytically in this work [52]. Equation (11) was solved following a LU (Lower -Upper triangular matrix) decomposition and back substitution afterwards. The above algorithm is enough, where the Jacobian matrix is well-conditioned. Otherwise, especially near the critical point, singular value decomposition was employed. With this provision, the algorithm switches smoothly from one incipient phase to the other by passing over the critical point.

RESULTS AND DISCUSSION
Compositions for four natural gas-like mixtures [53] are provided in Table 3. Figure 1 shows the phase envelopes calculated from PPR78, PR2SRK and NB for NG1 mixture. NG1 is a lean natural gas-like mixture with 91% methane. Compared to the experimental data, PPR78 markedly predicts the phase envelope while PR2SRK and NB are alike and a little inferior for dew point branch. The bubble point branch was equally well predicted by the three EoS. Figure 2 depicts the phase envelope for NG2 gas mixture. NG2 is a natural gas-like mixture with 22 components. NG2 contains 78% methane with the rest as branched and straight hydrocarbons with small amounts of carbon dioxide and aromatics. Shown in Figure 2, PPR78 is clearly in best agreement with the experimental data. PR2SRK and NB are alike, especially near the cricondentherm. Near the cricondenbar, these two a little separate but not appreciably. The phase envelopes for NG3 and NG4 are presented in the supplementary material. Figures S1 and S2 exhibit the phase envelopes, respectively. Similar results as mentioned above were found considering Figures S1 and S2. In general, PR2SRK and NB are alike and when compared to PPR78, they overpredict the cricondentherm and cricondenbar for natural gas-like mixtures. Real natural gas mixtures differ from the studied synthetic natural gas mixtures in a few ways. They may contain naphthenics and aromatics in addition to paraffins. Number of compounds comprising a real natural gas mixture may exceed 100. The composition of a real natural gas mixture is not often completely defined. Normally, few components up to C 5 or C 6 are determined and the rest is lumped and reported as unknown fraction, also called C 7þ . Unknown fractions are specified by two properties out of molecular weight, normal boiling point and specific gravity. Indeed, describing an unknown fraction as a component is oversimplification. It has been a common procedure to discretize a C 7þ by generalized SCN groups [54]. The more the number of SCN groups, the more accurate becomes equilibrium calculations. However, this costs computer runtime and memory. It was often suggested that for real natural gases and gas condensates, 12 SCN adequately describes  Phase envelope for natural gas-like mixture NG1 (experimental data from [53]). Phase envelope for natural gas-like mixture NG2 (experimental data from [53]).
unknown fractions [4]. The critical properties and acentric factor for SCN groups are determined using Twu's correlations [8] and Edmister formula [55], respectively. This work complies with this finding after employing exponential function provided by Pedersen et al. [56] to unknown fractions: Equation (13) is constrained to where z is the SCN composition, M is the SCN molecular weight and C and D are the two constants for the discretization model. C and D are determined for each unknown fraction by inserting Equation (13) into Equations (14), (15), omitting C and solving iteratively for D. Once D is determined, C can be calculated from Equation (14) without difficulty. Twelve real natural gas mixtures were examined in this work [57][58][59][60][61][62][63][64]. The compositions and C 7þ specifications are provided in Table S1 in the supplementary material. For predicting the binary interaction parameters, an effective structure for each SCN needed to be defined. SCN groups such as heptanes, octanes, and so on are complicated mixtures of paraffins, aromatics and naphthenics. In fact, SCN groups include all the components with normal boiling points within the range of two consecutive n-alkanes. A sophisticated method of calculating binary interaction parameters between a pair of compounds involving at least one SCN group was developed by Xu et al. [65]. Based on the molecular weight, normal boiling point or specific gravity of SCN groups, they determined the critical properties and acentric factor of the SCN groups. Then, assuming that the chemical structures of the SCN groups could merely be defined by paraphinic (CH 2 ), naphthenic (CH 2,cyclic ) and aromatic groups (CH aro ), Xu et al. [65] estimated the binary interaction parameters using Equation (4) and simultaneous solution of a system of three equations with three unknowns. This method is promising but a bit involved. A simpler model, however, can be proposed by assuming the structure of n-heptane for heptanes, n-octane for octanes and so on. Based on this idea, we performed the calculations with reasonable accuracy as shown afterwards. Figure 3 illustrates the phase envelope for RG1 mixture. RG1 is a sweet natural gas mixture. The experimental value of the dew point lies inside the phase envelopes and all three EoS overpredict the dew point value. The percent absolute deviation (%AD) in predicting the dew point pressure at the dew point temperature was found to be 14.3%, 23.8% and 20% using PPR78, PR2SRK and NB, respectively. We found no critical point. While the dew point curve extends from 200 K to 500 K and rises to very high pressure close to 200 K, the bubble point curve is short and extends from near 100 K to 200 K. The phase envelope shows a discontinuity near 200 K.
Exhibited in Figure 4a is the phase envelope for RG2 mixture. RG2 is a sour rich natural gas with 8.6% carbon dioxide. The dew point branches predicted by the three EoS underestimate the experimental value. The %AD in predicting the dew point pressure was found to be 20.2%, 14.8% and 17.9% using PPR78, PR2SRK and NB, respectively. Without employing binary interaction parameters, Nasrifar et al. [4] found a %AD of 30.36% in predicting the dew point pressure using PR78 and 26.18% using SRK. Clearly, the use of binary interaction parameters greatly improves the predictive capability of the EoS for highly sour natural gas mixtures and must not be avoided. The bubble point branch is peculiar, however. As can be seen the bubble point curve ends up at the critical point. For clarity the enlarged bubble point curve is shown in Figure 4b. Figure 4b demonstrates the presence of a second liquid phase l 2 . This was not unprecedented [6] knowing the high amount of carbon dioxide in RG2. Liquid rich carbon dioxide form a new phase. The presence of solid waxes cannot also be ruled out knowing that gas condensates are amenable to form waxes [66]. However, in this work it was not possible to predict solid wax equilibria. Figure 5 exhibits the two-phase envelope for RG3. The two EoS adequately predict the experimental value at dense phase region; however, NB with a %AD of 1.6% is slightly superior. The corresponding %AD for PPR78 was found to Phase envelope for real natural gas nature RG1 (experimental data from [60]). be 3.7%. Figure 6 is the same as Figure 5 while PPR78 and PR2SRK were used for calculating the phase envelope for RG4. Both EoS satisfactorily predict the phase envelope. The %AD was found to be 2.3% and 7.8% for PPR78 and PR2SRK, respectively. The phase envelopes for RG5 to RG12 are similar and provided in the supplementary material through Figures S3-S10. As reported by Nasrifar et al. [4], everything the same, the percent average absolute deviation (%AAD) in predicting the dew points for the same real natural gas mixtures reported in Table S1 without employing binary interaction parameters amounted to 13.4% for PR78 and 11.6% for SRK. While the error may look like high, it is a common practice to alleviate this inaccuracy by shifting the phase envelope towards the available experimental data and exploit the other EoS capabilities such as potential liquid hydrocarbon dropout in process design [5]. Presented in Table 4, the use of temperature-dependent binary interaction parameters clearly improved the PPR78 in predicting phase envelopes for gas mixtures with unknown fractions. The %AAD in predicting the dew point pressures is reported to be 11.2%; better than 13.4% documented in our previous work [4] employing zero binary interaction parameter and the same characterization method [8]. Indicated further in Table 4, among the EoS studied in this work, NB is slightly more accurate than the other two EoS. (a) Phase envelope for real natural gas nature RG2 (experimental data from [59]) and (b) Lower branches of the phase envelope for real natural gas nature RG2 (experimental data from [59]). Phase envelope for real natural gas mixture RG3 (experimental data from [64]). Phase envelope for real natural gas mixture RG4 (experimental data from [62]).

CONCLUSION
This work asserts that the use of temperature-dependent binary interaction parameters improved the applications of PPR78 in phase envelop calculations. Overall, the PR2SRK remained unchanged and NB got better results than the other EoS. The presence of a second liquid phase can be detected theoretically if a natural gas mixture contains high amount of carbon dioxide. This can be found by liquid-liquid equilibrium calculations.
Near the critical point, the Jacobian matrix becomes illconditioned. This can be used as a means to detect the critical point.  The %AAD in predicting dew points for natural gas mixtures reported in Table S1 using the studied EoS with and without binary interaction parameter.