Comparison of modifi ed alpha functions of the PR-EoS for volume and enthalpy prediction of natural gases

Proper simulation of processes of the natural gas industry such as dehydration, liquefaction and regasifi cation require accurate prediction of thermodynamic properties of the working fl uids. For such processes, cubic equations of state are the calculation methods most frequently employed. Among them, the Peng-Robinson equation is usually the one recommended for gas, refi nery and petrochemical applications in many simulators. Numerous works have been proposed in order to improve the temperature dependence relation of the attraction parameter of the equation – the so called alpha function. In this work, fi ve currently available alpha functions are evaluated for the prediction of molar volumes and enthalpies of natural gas samples. Additionally, parameters of one of the models are readjusted to volumetric data of methane, in order to represent its supercritical behavior more accurately. Experimental data of 44 mixtures are compared with calculated results. Van der Waals mixing rules are used, with binary interaction parameters set as zero. In the case of the original alpha function, it is also tested how the inclusion of non-zero binary parameters affects the predictions. The extended Saffari-Zahedi model presents the smallest average deviation for the molar volumes (1.35%). For the enthalpy calculation, the inclusion of the binary parameters results in deviation values of 2.62% for gas-gas transitions and 4.44% for gas-liquid transitions.


INTRODUCTION
Increasing global energy demand and growing concern with environmental issues accelerate the development of clean and economical energy sources. Low greenhouse gas and air pollutant emissions make the natural gas an attractive possibility, in comparison with other fossil fuels (Qyyum et al. 2018). The U.S Energy Information Administration expects the natural gas global consumption to grow 1.9% annually. Projections of the same agency indicate that by 2030 natural gas will surpass coal as the second most consumed fuel (U.S. EIA 2016).
For long-distance trading, natural gas is liquefi ed and stored in cryogenic tanks of ships (Wang et al. 2011). The design and optimization of liquefaction processes are usually conducted with software packages like Aspen Hysys, Aspen Plus and Honeywell UniSim Design, which employ accurate thermodynamic models for the estimation of physical properties of natural gas and refrigerant fl uids (Yuan et al. 2015). Saffari & Zahedi (2013) consider that the favorable balance between precision, simplicity and computational time of cubic equations of state (EoSs) -especially Peng-Robinson (PR), Soave-Redlich-Kwong (SRK) and Lee-Kesler-Plocker (LPL) -make them the calculation methods most frequently used in practical applications. According to the authors, the PR EoS is probably the most popular.
Improvements of thermodynamic properties predictions employing EoSs can be achieved by proper parameter adjustment for pure components or improvements in mixing rules, when mixtures are involved (Novak et al. 2018). For pure components a better performance of the equation can be obtained through refinements in the temperature dependence relation of the term of attraction -the so-called alpha function -by adjusting the function parameters to experimental data of vapor pressures, pressurevolume-temperature (PVT) data (Rijkers 1991, Trebble & Bishnoi 1987 or fugacity data (Flöter et al. 1998). Gasem et al. (2001) proposed modifications in the original PR alpha function and obtained 1.1% mean deviation in the vapor pressure prediction of 28 pure components, including heavy hydrocarbons (> C10), totaling 1100 experimental points. Saffari & Zahedi (2013) proposed a new alpha function for the Peng-Robinson equation, comparing calculated and experimental vapor pressure data of 10 components usually found in natural gases, from the triple point to the critical point. The new model presented a mean deviation of 1. 42%. Coquelet et al. (2004) compared the functions of Mathias & Copeman (1983), Trebble & Bishnoi (1987) and original PR (Peng & Robinson 1976), as well as proposing a new model, based on the first two equations. They evaluated 22 components and found an average absolute deviation of 1.2% for vapor pressure predictions, compared to 2.1% of the original PR equation.
Some researchers also reported how alpha functions affect phase equilibrium calculations for mixtures (Danesh et al. 1995, Twu et al. 1996, Flöter et al. 1998). However, works in which the models are evaluated for the estimation of volumetric or calorific properties of mixtures are scarce (Li et al. 2016). In order to study how alpha functions influence the prediction of these propertie, especially for natural gases, 5 different models are analyzed: PR's original alpha function, Gasem (Gasem et al. 2001), Mathias-Copeman (Mathias & Copeman 1983), Saffari-Zahedi (Saffari & Zahedi 2013) and Coquelet (Coquelet et al. 2004). The models are tested for the prediction of molar volumes of 10 natural gas samples and enthalpy variation of 34 mixtures of species commonly found in natural gases. Classical mixing rules are used, in which the binary interaction parameters (BIPs) are set equal to zero for all components (k ij = 0). The results are compared with the original PR EoS with non-zero BIPs.

THEORETICAL BACKGROUND: THERMODYNAMIC PROPERTIES DERIVED FROM EQUATIONS OF STATE
Equations of state together with models for ideal-gas heat capacities enable the calculation of physical properties of fluids (Vestfálová 2017). Using a modified version of Benedict-Webb-Rubin's equation, fundamental relations of classical thermodynamics and properly chosen thermodynamic paths, Younglove & Ely (1987) calculated the entropy, enthalpy, internal energy, heat capacities at constant pressure and volume, and speed of sound for methane, ethane, propane, isobutane and n-butane. A similar methodology was used by Vestfálová (2017) to calculate the thermal properties of water vapor. By applying mixing rules, one can extend the calculations for mixtures (Terron 2017).
The PR EoS (Peng & Robinson 1976) can be written as (equation 1): where P is the pressure, T is the temperature, v is the molar volume and R is the universal gas constant. Parameters a(T) and b result from mathematical restrictions on the critical point (equations 2 and 3): where subscript c indicates critical property and α(T) represents the alpha function.
When written in terms of compressibility factor (z), equation 1 becomes: The PR alpha function, originally proposed by Soave (1972) has the form (equations 5 and 6): where ω is the acentric factor. Parameter κ is generated by fitting vapor pressure data of a set of pure compounds. Equation 6 relates κ to the acentric factor of each substance. When mixtures are involved, mixing rules must be employed for the calculation of the parameters a(T) and b. In this work, the classical mixing rules of Van der Waals were adopted, as follows (equations 7, 8 and 9): The term k ij is the binary interaction parameter, specific for each pair of components of the mixture. They are usually adjusted from vapor-liquid equilibrium data, as described by Abudour et al. (2014).
The specific enthalpy variation (Δh) of a gas that undergoes a temperature and volume change from T 1 to T 2 and v 1 to v 2 is given by equation 10: where the abbreviation res relates to residual properties and Cp id is the heat capacity of the ideal gas at constant pressure. The specific residual Helmholtz free energy (a res ) offers a convenient possibility for deriving properties such as residual enthalpy and heat capacity at constant pressure, as described by Neau et al. (2009). It can be calculated as follows (equation 11): From the fundamentals of classic thermodynamics, it is possible to relate the residual Helmholtz free energy to specific residual enthalpy (equation 12): Subscript v denotes that the specific volume is kept constant during differentiation.
By coupling equations 11 and 12, an expression for the residual enthalpy as a function of compressibility factor is obtained (equation 13): Using equation 4 for the compressibility factor and solving the resulting integral, the specific residual enthalpy becomes (equation 14): Equation 15 results from rewriting the first integral of equation 10 in terms of the ideal-gas heat capacities of the components, weighted by their mole fractions. Once the residual enthalpies are computed using expression 14, the determination of the enthalpy variation for a mixture is straightforward (equation 15).

Alpha functions
The alpha functions used in molar volume and enthalpy calculations are shown in Table I. It must be noted that the models of Mathias & Copeman (1983) and Coquelet et al. (2004) make use of switching functions in order to represent supercritical behavior more accurately. The other 3 models consist of one function each. The terms Ci, ci and k i are constants for each component. A, B and C are constants of the Gasem model and T r is the reduced temperature.
The Saffari-Zahedi model was originally proposed under saturation conditions. In order to evaluate how a better adjustment of the alpha function for methane, only, would impact the results, the parameters of the same function were readjusted to the data generated by Setzmann & Wagner (1991), based on a fundamental equation explicit in Helmholtz free energy containing 40 fitted coefficients. The authors reported deviations from experimental densities of 0.03% for pressures below 12 MPa and temperatures below 350K and deviations of 0.03% to 0.15% for higher pressures and temperatures. A total of 200 points in the temperature range of 195 K to 430 K and pressure range of 200 kPa to 8000 kPa was used. The first parameter k 1 was kept the same value as in the original function. By doing so, it is possible to avoid discontinuities in the alpha function. The two remaining parameters, k 2 and k 3 were recalculated by minimizing the objective function f -sum of the square of relative density (ρ) deviations -, given by equation 16: The subscripts calc and exp relate to calculated and experimental values, respectively. Parameter N exp is the number of experimental points used.
The optimization code was developed in the software MATLAB. The function fminsearch was used. Parameters of the original and extended models for methane are given in Table II, where T r is the reduced temperature.

Data source
PVT data of 10 natural gases (GN) were obtained from different sources in the literature. Details of the data sources, compositions and ranges of temperature and pressure are given in Tables  III and IV. Uncertainties of volume and density measurements were only found for mixture GN10, for which uncertainties of 0.044% are estimated for the first 4 experimental points (Magee et al. 1997). For the remaining points (Simon et al. 1997), uncertainties of 0.04% are estimated.
Experimental enthalpy data were obtained from Ashton & Haselden (1980). The authors have measured the enthalpy variation of 34 mixtures containing different mole fractions of nitrogen, methane, ethane, propane, isobutane and n-butane, undergoing a cooling process at constant pressure. A total of 8 points includes transitions between gas and sub-cooled liquid regions. The total uncertainty of enthalpy

Function
Reference/Author Once z was computed, the molar volumes were calculated with equation 20:  Enthalpy variations were obtained from equations 14 and 15.
When applying the mixing rules, the binary interaction parameters were set equal to zero (k ij =0). For the PR original alpha function, it was also evaluated how the inclusion of BIPs would affect the predictions. In that case, BIPs were obtained from Li et al. (2016) and considered independent of temperature. The binary parameters for helium were considered zero since no data were available. Critical pressures and temperatures of the components and correlations and parameters for calculating ideal-gas heat capacities were taken from Terron (2017). The prediction accuracy for both molar volume and enthalpy was measured where N exp is the number of experimental points.

RESULTS
The volumes of the natural gas samples -GN1 to GN10 -were calculated by solving equation 1 for v. Table VI comprises calculated volumes for the mixture GN10, as an example. The smallest deviations were obtained with the Saffari-Zahedi extended model, as seen in Table VI. Apart from that model, for the first four experimental points, the inclusion of the BIPs resulted in a small improvement in the prediction of specific volumes. The effects associated with the difference in size and polarity between the molecules of hydrocarbons and the molecules of nitrogen and carbon dioxide seem to become more intense as the temperature decreases. Up from 249.996 K the effect of including the BIPs The mixtures GN3, GN4, GN5 and GN6 differ essentially by the concentration of carbon dioxide, being all under the same temperature, 310 K. As this concentration increases, an increase in the deviations of the models for which the binary interaction parameters are equal to zero becomes evident. This trend was also observed by Li et al. (2016), who suggest that large deviations between experimental and calculated compressibility factors may occur when the PR EoS is applied to high CO 2 -content natural gas mixtures, especially when the CO 2 mole fraction exceeds 10%. Only for the mixture GN7, whose concentration of carbon dioxide exceeds 50% and methane concentration is less than 45%, the effect of including BIPs becomes more significant than the adjustment of the alpha function proposed for methane, as seen in figures 1a, 1b, 1c and 1d. Figures 1a and 1b illustrate the behavior of the original PR EoS and the Coquelet model, respectively. Figure 1c and 1d represent the PR EoS with non-zero BIPs and extended Saffari-Zahedi model, in that order.
The effects of a temperature increase can be analyzed with the data of the mixtures GN1, GN8 and GN9. For GN1 an increase in temperature results in reduction of deviations for 4 of the 5 original models for each pressure value analyzed (figures 2a, 2b, 2c and 2d). The Gasem model do not follow a similar trend (figure 2e). For the extended model, there is a pressure range (100kPa-250kPa) in which an increase in temperature causes the deviation to increase. Above this range, the trend is the same as that described for the first 4 models. The point within the curve of 444.26 K with an abnormally small deviation is believed to be an outlier.
In the case of the pair GN8-GN9, the same reduction trend can be verified for pressures  up from 50 bar ( figure 3). An exception is the extended model (figure 3f), for which deviations are slightly bigger for higher temperatures in the pressure range 0-120 kPa. Though the point corresponding to the pressure of 34.83 bar is expected to be an outlier, more experimental points in the pressure range from 20 bar to 50 bar would be needed in order to draw a conclusion. For enthalpy calculation, the same alpha functions were used. The models were employed for the calculation of the parameter a(T), in equation 14. Mixing rules given by equations from 7 to 9 were also used. The results from tables VI to X indicate that the inclusion of binary parameters in the mixing rules must be analyzed carefully, since the prediction accuracy is affected by temperature and concentration of incondensable components. When such components -in this work, N 2 and CO 2 -are present in high concentrations, the incorporation of BIPs results in improvements in the prediction of volumes and enthalpies. When not, adjusting the alpha function to volumetric data of supercritical methane produces better volume estimations.
Changing the alpha function, however, does not improve significantly the enthalpy predictions. Actually, adjusting the alpha function to density data of methane can even produce worst enthalpy predictions, as seen for the Saffari-Zahedi extended model. A possible explanation would be a poor representation of the derivatives of the alpha functions with temperature, despite a good adjustment of the alpha function itself, since the residual enthalpy depends on both the alpha function and on its first derivative, as seen in equation 14. However, an examination of expected values for the alpha function and for its derivative suggests this is  Figure 4a was generated by isolating the term corresponding to the alpha function in the PR EoS and calculating its value in order to represent the PVT data from Setzmann & Wagner (1991). In figure 4b, the derivatives for each pressure were calculated by using finite differences approach, with a step size of 10 K. Although the extended model is relatively well adjusted to the PVT-based curves (figure 4a) and the derivative presents the same tendency of monotonic increment as the other curves (figure 4b), using this model and its derivative led to worse predictions of enthalpy for the natural gas samples.
The inclusion of binary parameters has also a small effect on enthalpy calculation for processes in the gas region. When phase transitions are presented, however, their inclusion results in more meaningful improvements, as seen in Table X.
In comparison with the calculations performed by Ashton & Haselden (1980) and except for the Saffari-Zahedi extended model, the enthalpy average deviations related to gasgas transitions in this work were smaller than those calculated with the BWR (Benedict-Webb-Rubin) equation, with original and modified mixing rules -which resulted in deviations of 3.16% and 3.12%. Employing the generalized BWRS (Benedict-Webb-Rubin-Starling) equation with non-zero BIPs, the authors obtained an average deviation of 2.52%, which is slightly smaller than the average deviation when the PR EoS is applied with non-zero BIPs (2.62%). For gasliquid transitions, the biggest average deviation in this work -5.95% -was significantly smaller than that of the article previously mentioned -9%, in the best case. As the best result, an average deviation of 4.44% was obtained with the PR-BIPs model.

CONCLUSIONS
Simulation and design of processes commonly found in the natural gas industry require accurate prediction of thermodynamic properties. For this purpose, cubic equations of state are commonly used in the industry, given the favorable balance between precision and computational time they offer. In this work, it was investigated if a better performance of the Peng-Robinson equation could be achieved by improving the temperature dependence relation -or alpha function -of its attraction parameter.
Based on experimental data, five different alpha functions were evaluated for the prediction of molar volume and enthalpy of natural gas samples. For methane, the Saffari-Zahedi model   originally proposed under saturation conditions was readjusted to volumetric data above the critical temperature. The readjusted function presents the best performance for estimation of molar volumes, with an average AAD of 1.35%. AADs of 2.62% for enthalpy variations of transitions within the gas region and 4.44% for enthalpy variations of gas-liquid transitions resulted when the original Peng-Robinson equation was employed with non-zero binary parameters (BIPs), whose inclusion proved to be relevant as the concentration of incondensable components -N 2 and CO 2 -increased. Improving the alpha function has not led to significant improvements in enthalpy calculations for the natural gas samples. Although the strategy of adjusting the alpha function to volumetric data of methane resulted in significantly better density predictions for the mixtures -producing results even better than those obtained by just including binary parameters in the original equation -, it has worsened the enthalpy predictions.