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

Abstract

Proper simulation of processes of the natural gas industry such as dehydration, liquefaction and regasification require accurate prediction of thermodynamic properties of the working fluids. 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, refinery 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, five 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.

Key words
Alpha function; enthalpy; equation of state; natural gas; Peng-Robinson; volume

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 liquefied 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 fluids (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, pressure-volume-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):

$$P=\frac{RT}{v-b}-\frac{a(T)}{v(v+b)+b(v-b)}$$ (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):

$$a(T)=\left[\frac{0.45724R^2{T}_c^2}{P_c^2}\right]\alpha (T)$$ (2)

$$b=\frac{0.0780R{T}_c}{P_c}$$ (3)

where subscript c indicates critical property and α(T) represents the alpha function.

When written in terms of compressibility factor (z), equation 1 becomes:

$$z=\frac{v}{v-b}+\frac{-a(T)/RT}{v(v+b)+b(v-b)}$$ (4)

The PR alpha function, originally proposed by Soave (1972) has the form (equations 5 and 6):

$$a(T)={[1+\kappa (1-{\left(\frac{T}{T_c}\right)}^{0.5})]}^2$$ (5)

$$\kappa =0.37464+1.54226\omega -0.26992{\omega }^2$$ (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):

$$a(T)=\sum _i\sum _j{y}_i{y}_j{a}_{ij}$$ (7)

$$b=\sum _i^{}{x}_i{b}_i$$ (8)

where

$$a_{ij}=\sqrt{a(T{)}_{i}a(T{)}_j}(1-k_{ij})$$ (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₁ to T₂ and v₁ to v₂ is given by equation 10:

$$\mathrm{\Delta }{h}_{1\to 2}=\underset{T_1}{\overset{T_2}{\int }}C{p}^{id}dT+h_2^{res}(T_2,v_2)-h_1^{res}(T_1,v_1)$$ (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):

$$a^{res}=-RT{\int }_{\mathrm{\infty }}^v(z-1)\frac{dv}{v}$$ (11)

From the fundamentals of classic thermodynamics, it is possible to relate the residual Helmholtz free energy to specific residual enthalpy (equation 12):

$$h^{res}=RT(z-1)-T^2\left(\frac{\mathrm{\partial }{a}^{res}/T}{\mathrm{\partial }T}\right)$$ (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):

$$h^{res}=RT(z-1)+R{T}^2{\int }_{\mathrm{\infty }}^v\left({\frac{\mathrm{\partial }z}{\mathrm{\partial }T}}_v\right)\frac{dv}{v}$$ (13)

Using equation 4 for the compressibility factor and solving the resulting integral, the specific residual enthalpy becomes (equation 14):

$$h^{\mathit{res}}=\mathit{RT}\left(z-1\right)+\left[a\left(T\right)-T\frac{\mathit{da}\left(T\right)}{\mathit{dT}}\right]\mathrm{.}\left[\frac{1}{2\sqrt{2}b}\ln \left(\frac{v-b\left(\sqrt{2}-1\right)}{v+b\left(\sqrt{2}+1\right)}\right)\right]$$ (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).

$${∆h}_{1\to 2}=\underset{T_1}{\overset{T_2}{\int }}\sum _i{{\mathit{Cp}}^{\mathit{id}}}_i{y}_i\mathit{dT}+h_2^{\mathit{res}}(T_2,v_2)-h_1^{\mathit{res}}(T_1,v_1)$$ (15)

MATERIALS AND METHODS

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.

Table I
Alpha functions.

Function	Reference/Author
$$\alpha ={\left[1+S\left(1-\sqrt{T_r}\right)\right]}^2$$ $$S=0.37464+1.5422\omega -0.26992{\omega }^2$$	(Peng & Robinson 1976)
$$\alpha ={[1+C_1\left(1-\sqrt{T_r}\right)+C_2{\left(1-\sqrt{T_r}\right)}^2+C_2{\left(1-\sqrt{T_r}\right)}^3]}^2$$ $$C_2=C_3=0,T_r>1$$	(Mathias & Copeman 1983)
$$\alpha =\exp [\left(A+B{T}_r\right)\left(1-{T_r}^S\right)]$$ $$S=C+\mathit{D\omega }+E{\omega }^2$$	(Gasem et al. 2001)
$$\alpha =\exp \left[c_1\left(1-T_r\right)\right]{[1+c_2{\left(1-T_r^{1/2}\right)}^2+c_3{\left(1-T_r^{1/2}\right)}^3]}^2,T_r<1;$$ $$c_2=c_3=0,T_r>1$$	(Coquelet et al. 2004)
$$\alpha =\exp [k_1{T}_r+k_2\ln T_r+k_3\left(1-T_r^{1/2}\right)]$$	(Saffari & Zahedi 2013)

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₁ 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₂ and k₃ were recalculated by minimizing the objective function f – sum of the square of relative density (ρ) deviations –, given by equation 16:

$$f=\sum _i^{N_{\exp }}{\left(\frac{{\rho }_{\mathit{calc}}(i)-{\rho }_{\exp }(i)}{{\rho }_{\exp }(i)}\right)}^2$$ (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.

Table II
Parameters of the original and extended Saffari-Zahedi alpha function for methane.

$$\alpha =\exp [k_1{T}_r+k_2\ln T_r+k_3\left(1-T_r^{1/2}\right)]$$
k1	k2	k3	Tr
0.0039	0.0473	0.8514	≤1
0.0039	-0.1715	0.8979	>1

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.

Table III
Compositions of the mixtures used for molar volume prediction.

Components	Composition (mole fraction)
	GN1	GN2	GN3	GN4	GN5	GN6	GN7	GN8	GN9	GN10
CH4	0.8074	0.8272	0.9740	0.9222	0.8616	0.7830	0.4460	0.7559	0.6993	0.81299
C2H6	0.0870	0.1470	0.0153	0.0149	0.0146	0.0139	0.0068	0.0140	0.0106	0.03294
C3H8	0.0290	0.0104	0.0008	0.0009	0.0004	0.0004	-	0.0021	0.0005	0.00637
i-C4H10	-	0.0004	-	-	-	-	-	-	-	0.00101
n-C4H10	0.0170	0.0008	-	-	-	-	-	-	-	0.00100
n-C5H12	0.0013	-	-	-	-	-	-	-	-	-
n-C6H14	0.0006	0.0001	-	-	-	-	-	-	-	-
n-C7H16	0.0003	-	-	-	-	-	-	-	-	-
N2	0.0480	0.0062	0.0095	0.0084	0.0088	0.0055	0.0026	0.0050	0.0082	0.13575
CO2	0.0090	0.0079	0.0004	0.0536	0.1146	0.1972	0.5446	0.2230	0.2814	0.00994
He	0.0003	-	-	-	-	-	-	-	-	-
Number of points	52	10	5	5	6	6	8	6	6	23

Table IV
Reference, pressure and temperature range of the mixtures used for molar volume prediction.

Mixture	Pressure range (MPa)	Temperature range (K)	Reference
GN1	1.38-55.2	344.26-444.26	(Lee et al. 1966) apud (Terron 2017)
GN2	1.00-15.0	301.15	(Hannisdal 1987) apud (Terron 2017)
GN3	2.43-10.7	310.93	(Robinson Jr. & Jacoby 1965) apud (Terron 2017)
GN4	2.36-10.3	310.93
GN5	1.77-11.1	310.93
GN6	1.75-10.8	310.93
GN7	1.06-11.7	310.93
GN8	1.72-11.2	338.71
GN9	1.61-10.8	366.48
GN10	1.99-19.7	225.01-324.99	(Magee et al. 1997); (Simon et al. 1997)

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 measurements is reported to be 0.90%. Details can be found in Table V.

Table V
Composition, pressure and temperature range of the mixtures used for enthalpy variation prediction.

Mixtures	Pressure, bar	Temperature, K		Mole fraction						Final state
		Tin	Tout	N2	CH4	C2H6	C3H8	i-C4H10	n-C4H10
1	13.89	366.91	200.32	0.0445	0.9277	0.0231	0.0040	0.0003	0.0004	Gas
2	27.58	367.02	231.46	0.0510	0.9210	0.0233	0.0040	0.0003	0.0004
3	27.99	367.06	200.25	0.0522	0.9198	0.0233	0.0040	0.0003	0.0004
4	44.63	367.47	255.79	0.0562	0.9161	0.0232	0.0038	0.0003	0.0004
5	41.44	367.29	224.02	0.0575	0.9148	0.0232	0.0038	0.0003	0.0004
6	14.07	367.02	240.23	0.0739	0.8972	0.0232	0.0048	0.0004	0.0005
7	13.73	366.88	202.80	0.0752	0.8959	0.0232	0.0048	0.0004	0.0005
8	27.44	366.66	226.93	0.0768	0.8946	0.0232	0.0047	0.0004	0.0004
9	15.17	361.00	241.64	0.2685	0.7037	0.0201	0.0048	0.0014	0.0015
10	15.02	360.94	217.11	0.2697	0.7029	0.0197	0.0048	0.0014	0.0015
11	27.50	360.84	237.80	0.2736	0.7004	0.0186	0.0046	0.0013	0.0015
12	27.05	360.93	225.19	0.2740	0.7002	0.0185	0.0044	0.0013	0.0015
13	27.06	360.91	207.97	0.2752	0.6990	0.0185	0.0045	0.0013	0.0015
14	41.45	361.85	227.37	0.2716	0.7001	0.0207	0.0047	0.0014	0.0015
15	42.16	361.80	213.20	0.2720	0.7001	0.0205	0.0045	0.0014	0.0015
16	41.62	361.85	207.78	0.2745	0.6976	0.0205	0.0046	0.0014	0.0015
17	14.53	361.89	220.12	0.2005	0.7672	0.0254	0.0056	0.0006	0.0007
18	27.90	361.72	204.65	0.2107	0.7575	0.0248	0.0054	0.0007	0.0009
19	42.80	361.69	208.03	0.2096	0.7586	0.0257	0.0053	0.0004	0.0004
20	55.07	361.71	215.62	0.2142	0.7489	0.0277	0.0074	0.0010	0.0010
21	55.12	362.02	217.26	0.2147	0.7486	0.0272	0.0077	0.0010	0.0010
22	14.39	362.31	224.90	0.1576	0.8030	0.0302	0.0072	0.0009	0.0011
23	14.55	362.28	205.60	0.1580	0.8026	0.0302	0.0072	0.0009	0.0011
24	27.48	361.67	215.84	0.1585	0.8017	0.0304	0.0074	0.0010	0.0011
25	41.09	361.32	215.21	0.1527	0.8076	0.0305	0.0070	0.0011	0.0011
26	55.66	361.44	212.64	0.1532	0.8074	0.0305	0.0067	0.0011	0.0011
27	41.39	367.00	128.86	0.0594	0.9130	0.0231	0.0038	0.0003	0.0004	Sub-cooled liquid
28	41.08	366.91	167.80	0.0602	0.9122	0.0231	0.0038	0.0003	0.0004
29	41.25	361.91	161.97	0.2739	0.6982	0.0205	0.0045	0.0014	0.0015
30	55.01	361.09	179.54	0.2691	0.6980	0.0241	0.0070	0.0009	0.0009
31	55.21	361.68	177.62	0.2146	0.7488	0.0273	0.0075	0.0010	0.0010
32	42.02	361.12	173.71	0.1525	0.8077	0.0306	0.0070	0.0011	0.0011
33	56.46	360.92	189.55	0.1535	0.8071	0.0305	0.0067	0.0011	0.0011
34	56.27	361.01	184.41	0.1532	0.8074	0.0305	0.0067	0.0011	0.0011

Calculation and accuracy analysis

Enthalpy and molar volumes calculations were conducted in the software MATLAB. The compressibility factors (z) were computed by rewriting the PR EoS in a polynomial form (equations 17, 18 and 19). The function vpasolve was used for solving equation 17. Parameters a(T) and b were calculated with Van der Waals mixture rules (equations 7, 8 and 9).

$$z^3-{\left(1+B-2B\right)z}^2+(A-2B-3B^2)z+(\mathit{AB}-B^2-B^3)=0$$ (17)

$$A=\frac{a(T)P}{{(\mathit{RT})}^2}$$ (18)

$$B=\frac{\mathit{bP}}{\mathit{RT}}$$ (19)

Once z was computed, the molar volumes were calculated with equation 20:

$$v=\frac{\mathit{RTz}}{P}$$ (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 by calculating the average absolute deviation (%AAD) - given by equations 21 and 22 – and maximum absolute deviation (max. %AD) – given by equations 23 and 24.

$$\text{\%}\mathit{AAD}=\frac{1}{N_{\exp }}\times \sum _i^{N_{\exp }}\left(\frac{\left|{∆h}_{\mathit{calc}}(i)-{∆h}_{\exp }(i)\right|}{{∆h}_{\exp }(i)}\right)\times 100\text{\%}$$ (21)

$$\text{\%}\mathit{AAD}=\frac{1}{N_{\exp }}\times \sum _i^{N_{\exp }}\left(\frac{\left|v_{\mathit{calc}}(i)-v_{\exp }(i)\right|}{v_{\exp }(i)}\right)\times 100\text{\%}$$ (22)

$$\mathit{max}\mathrm{.}\text{\%}\mathit{AD}=\underset{i}{\mathit{max}}\left(\frac{\left|{∆h}_{\mathit{calc}}(i)-{∆h}_{\exp }(i)\right|}{{∆h}_{\exp }(i)}\right)\times 100\text{\%}$$ (23)

$$\mathit{max}\mathrm{.}\text{\%}\mathit{AD}=\underset{i}{\mathit{max}}\left(\frac{\left|v_{\mathit{calc}}(i)-v_{\exp }(i)\right|}{v_{\exp }(i)}\right)\times 100\text{\%}$$ (24)

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.

Table VI
Experimental and calculated molar volumes for mixture GN10.

Experimental data				Calculated volumes/ m³/kmol
				Equation of state/Alpha function
T/K	P/KPa	v/m³/kmol		PR-original	PR-Coquelet	PR-Saffari-Zahedi	PR-BIPs	PR-Saffari-Zahedi ext.
225.009	3193	0.4988		0.4865	0.4879	0.4850	0.4880	0.4901
225.007	4512	0.3248		0.3141	0.3156	0.3123	0.3158	0.3181
225.008	14818	0.07124		0.06911	0.06990	0.06827	0.06998	0.07118
225.009	19712	0.06213		0.05894	0.05942	0.05844	0.05946	0.06019
249.996	10317	0.14069		0.1365	0.1387	0.1346	0.1379	0.1418
250.009	8290	0.1864		0.1802	0.1824	0.1781	0.1816	0.1857
250.01	5967	0.2823		0.2736	0.2757	0.2716	0.2749	0.2788
250.008	3746	0.4889		0.4779	0.4798	0.4762	0.4791	0.4827
250.009	1986	0.9814		0.9688	0.9705	0.9672	0.9699	0.9732
274.994	10407	0.17513		0.1701	0.1724	0.1680	0.1710	0.1753
274.991	8654	0.2172		0.2109	0.2133	0.2088	0.2119	0.2163
274.992	6042	0.3291		0.3206	0.3230	0.3185	0.3216	0.3260
274.994	4167	0.4985		0.4883	0.4906	0.4863	0.4893	0.4935
274.994	2351	0.9217		0.9102	0.9124	0.9083	0.9111	0.9152
299.996	8569	0.2560		0.2493	0.2519	0.2471	0.2501	0.2545
299.998	5866	0.3877		0.3793	0.3820	0.3771	0.3801	0.3846
299.999	3984	0.5872		0.5778	0.5803	0.5755	0.5785	0.5829
299.996	2251	1.068		1.058	1.060	1.055	1.058	1.063
324.986	8812	0.2804		0.2741	0.2769	0.2717	0.2748	0.2790
324.988	5961	0.4246		0.4170	0.4198	0.4146	0.4176	0.4220
324.99	4820	0.5313		0.5229	0.5257	0.5204	0.5235	0.5279
324.986	4014	0.6435		0.6344	0.6373	0.6320	0.6351	0.6395
324.986	2162	1.220		1.209	1.212	1.207	1.210	1.2140
Average %AAD:				2.339	1.653	2.964	1.976	0.8778

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 is surpassed by the better adjustment of the Coquelet model.

The average deviations between the experimental and predicted volumes for all mixtures and the alpha function used in the predictions are shown in Table VII.

Table VII
Deviations of experimental and calculated molar volumes.

Mixture	Equation of state/Alpha function
	PR-original (%AAD)	PR-Mathias-Copeman (%AAD)	PR-Gasem (%AAD)	PR-Saffari-Zahedi (%AAD)	PR-Coquelet (%AAD)	PR-BIPs (%AAD)	PR-Saffari-Zahedi ext. (%AAD)
GN1	7.063	6.651	9.481	7.942	6.011	6.918	5.378
GN2	1.903	1.483	1.695	2.591	1.040	1.769	1.054
GN3	1.680	1.366	1.517	2.169	1.023	1.666	0.3827
GN4	1.976	1.675	1.820	2.457	1.346	1.677	0.4140
GN5	2.138	1.854	1.991	2.609	1.543	1.560	0.5380
GN6	2.329	2.055	2.188	2.804	1.755	1.404	0.7857
GN7	3.528	3.303	3.412	4.088	3.054	1.645	2.338
GN8	1.939	1.709	1.796	2.416	1.412	1.197	0.7145
GN9	1.694	1.558	1.593	2.054	1.352	1.218	1.025
GN10	2.339	2.015	7.090	2.964	1.653	1.976	0.8778
Average %AAD:	2.659	2.367	3.258	3.209	2.019	2.103	1.351

Table VIII presents the maximum absolute deviation of each mixture, employing different alpha-functions.

Table VIII
Maximum deviations between experimental and calculated molar volumes.

Mixture	Equation of state/Alpha function
	PR-original (max. %AD)	PR-Mathias-Copeman (max. %AD)	PR-Gasem (max. %AD)	PR-Saffari-Zahedi (max. %AD)	PR-Coquelet (max. %AD)	PR-BIPs (max. %AD)	PR-Saffari-Zahedi ext. (%AAD)
GN1	12.98	12.56	17.59	13.78	11.97	12.79	10.95
GN2	2.661	2.183	2.399	3.684	1.729	2.477	3.085
GN3	2.289	1.709	1.985	3.207	1.264	2.263	1.008
GN4	2.848	2.278	2.553	3.762	1.681	2.282	0.8979
GN5	3.470	2.855	3.151	4.493	2.183	2.221	0.7953
GN6	4.316	3.713	4.005	5.369	3.055	2.286	1.236
GN7	9.049	8.793	8.917	9.684	8.509	6.889	7.688
GN8	3.466	2.970	3.158	4.500	2.330	1.869	1.039
GN9	2.744	2.369	2.464	3.740	2.313	2.432	2.602
GN10	5.125	4.724	14.73	5.935	4.357	4.286	3.121
Average max. %AD:	4.895	4.415	6.096	5.815	3.934	3.979	3.246

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₂-content natural gas mixtures, especially when the CO₂ 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.

Figure 1
Volume deviations (%AAD) calculated with the models PR (a), Coquelet (b), PR with non-zero BIPs (c) and extended Saffari-Zahedi (d).

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.

Figure 2
Volume deviations (%AAD) of mixture GN1 calculated with the models PR with non-zero BIPs (a), Coquelet (b), Saffari-Zahedi (c), Mathias-Copeman (d), Gasem (e) and extended Saffari-Zahedi (f).

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.

Figure 3
Volume deviations (%AAD) of mixtures GN8 and GN9, calculated with the models PR with non-zero BIPs (a), Coquelet (b), Saffari-Zahedi (c), Mathias-Copeman (d), Gasem (e) and extended Saffari-Zahedi (f).

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. Table IX comprises the deviations for mixtures with no phase change while Table X contains results only for gas-liquid transitions.

Table IX
Deviations of experimental and calculated enthalpy variation for gas-gas transitions.

Mixture	Equation of state/Alpha function
	PR-original (%AD)	PR-Mathias-Copeman (%AD)	PR-Gasem (%AD)	PR-Saffari-Zahedi (%AD)	PR-Coquelet (%AD)	PR-BIPs (%AD)	PR-SZ ext. (%AD)
1	0.9089	1.014	0.9028	0.9018	1.008	0.8919	1.803
2	2.712	2.860	2.705	2.671	2.933	2.676	3.838
3	0.7339	0.9474	0.7275	0.7312	0.9402	0.677	2.535
4	5.434	5.608	5.434	5.361	5.778	5.374	6.636
5	1.679	1.898	1.671	1.651	1.986	1.598	3.378
6	3.350	3.424	3.348	3.321	3.470	3.331	3.885
7	0.2213	0.1225	0.2256	0.2322	0.1209	0.2467	0.5875
8	7.297	7.456	7.292	7.256	7.528	7.242	8.503
9	3.185	3.258	3.197	3.141	3.324	3.125	3.613
10	1.984	2.073	1.992	1.949	2.117	1.913	2.551
11	4.937	5.068	4.957	4.866	5.180	4.807	5.722
12	3.078	3.221	3.095	3.018	3.312	2.938	3.974
13	1.543	1.712	1.558	1.499	1.768	1.375	2.673
14	6.393	6.598	6.419	6.325	6.731	6.137	7.696
15	4.083	4.321	4.104	4.045	4.409	3.764	5.701
16	5.027	5.282	5.049	5.003	5.350	4.677	6.797
17	4.074	4.160	4.078	4.043	4.201	4.019	4.676
18	0.7999	0.9842	0.8079	0.7708	1.020	0.6420	2.149
19	1.416	1.661	1.427	1.414	1.724	1.104	3.351
20	0.6351	0.9105	0.6463	0.6777	0.9813	0.1815	2.657
21	1.176	1.448	1.188	1.208	1.528	0.7385	3.144
22	3.491	3.578	3.492	3.462	3.614	3.446	4.084
23	1.209	1.314	1.210	1.191	1.325	1.156	1.964
24	3.067	3.244	3.071	3.033	3.290	2.951	4.302
25	2.052	2.302	2.057	2.041	2.356	1.831	3.835
26	1.6034	1.912	1.602	1.713	1.934	1.151	3.926
%AAD:	2.773	2.938	2.779	2.751	2.997	2.615	3.954

Table X
Deviations of experimental and calculated enthalpy variation for gas-liquid transitions.

Mixture	Equation of state/Alpha function
	PR-original (%AD)	PR-Mathias-Copeman (%AD)	PR-Gasem (%AD)	PR-Saffari-Zahedi (%AD)	PR-Coquelet (%AD)	PR-BIPs (%AD)	PR-SZ ext. (%AD)
27	1.604	0.9979	1.200	1.049	1.551	1.797	1.289
28	0.3467	0.1941	0.1988	0.2499	0.06074	0.6240	0.1232
29	7.236	7.853	7.437	7.350	7.722	5.949	7.396
30	5.779	6.414	5.886	6.336	6.159	3.289	6.722
31	3.408	4.106	3.553	3.587	3.968	2.077	3.902
32	9.381	10.13	9.558	9.518	9.957	8.441	9.711
33	8.728	9.511	8.788	9.434	9.131	6.796	10.04
34	7.733	8.552	7.854	7.991	8.343	6.522	8.442
%AAD:	5.527	5.970	5.559	5.689	5.862	4.437	5.953

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₂ and CO₂ – 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 not a good hypothesis. Both are represented in figure 4. 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.

Figure 4
PVT-based and calculated values for the Peng-Robinson (PR), Safarri-Zahedi (SZ) and extended Saffari-Zahedi (SZext.) alpha functions (a) and its derivatives with temperature (b) for different pressures.

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 gas-gas 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 gas-liquid 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₂ and CO₂ - 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.

ACKNOWLEGMENTS

The authors are grateful to Conselho Nacional do Desenvolvimento Científico e Tecnológico (CNPq) – and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG, TEC-APQ-00914-16) – for the financial support.

REFERENCES

• ABUDOUR AM, MOHAMMAD SA, ROBINSON RL & GASEM KA. 2014. Generalized binary interaction parameters for the Peng–Robinson equation of state. Fluid Phase Equilibr 383: 156-173.
• ASHTON GJ & HASELDEN GG. 1980. Measurements of enthalpy and phase equilibrium for simulated natural gas mixtures and correlation of the results by a modified Starling equation. Cryogenics 20: 41-47.
• COQUELET C, CHAPOY A & RICHON D. 2004. Development of a new alpha function for the Peng-Robinson equation of state: Comparative study of alpha function models for pure gases (natural gas components) and water-gas systems. Int J Thermphys 25: 133-158.
• DANESH A, XU DH, TEHRANI DH & TODD AC. 1995. Improving predictions of equation of state by modifying its parameters for super critical components of hydrocarbon reservoir fluids. Fluid Phase Equilibr 112: 45-61.
• FLÖTER E, DE LOOS TW & DE SWAAN ARONS J. 1998. Improved modeling of the phase behavior of asymmetric hydrocarbon mixtures with the Peng-Robinson equation of state using a different temperature dependency of the parameter. Ind Eng Chem Res 37: 1651-1662.
• GASEM KA, GAO W, PAN Z & ROBINSON RL. 2001. A modified temperature dependence for the Peng-Robinson equation of state. Fluid Phase Equilibr 181: 113-125.
• HANNISDAL NE. 1987. Gas Compression Equations Evaluated. Oil Gas J: 38-41.
• LEE AL, GONZALEZ MH & EAKIN BE. 1966. The Viscosity of Natural Gases. J Petrol Technol 18: 997-1000.
• LI Z, JIA W & LI C. 2016. An improved PR equation of state for CO2-containing gas compressibility factor calculation. J Nat Gas Sci Eng 36: 586-596.
• MAGEE JW, HAYNES WM & HIZA MJ. 1997. Isochoric (p,ρ,T) measurements for five natural gas mixtures from T = (225 to 350) K at pressures to 35 MPa. J Chem Thermodyn 29: 1439-1454.
• MATHIAS PM & COPEMAN TW. 1983. Extension of the Peng-Robinson equation of state to complex mixtures: Evaluation of the various forms of the local composition concept. Fluid Phase Equilibr 13: 91-108.
• NEAU E, RASPO I, ESCANDELL J, NICOLAS C & HERNÁNDEZ-GARDUZA O. 2009. The Soave, Twu and Boston-Mathias alpha functions in cubic equations of state. Part II. Modeling of thermodynamic properties of pure compounds. Fluid Phase Equilibr 276: 156-164.
• NOVAK N, LOULI V, SKOURAS S & VOUTSAS E. 2018. Prediction of dew points and liquid dropouts of gas condensate mixtures. Fluid Phase Equilibr 457: 62-73.
• PENG DY & ROBINSON DB. 1976. A New Two-Constant Equation of State. Ind Eng Chem Fundam 15: 59-64.
• QYYUM MA, QADEER K & LEE M. 2018. Comprehensive Review of the Design Optimization of Natural Gas Liquefaction Processes: Current Status and Perspectives. Ind Eng Chem Res 57: 5819-5844.
• RIJKERS RP. 1991. Retrograde condensation of lean natural gas. Doctoral dissertation, Delft University of Technology. Available on: http://resolver.tudelft.nl/uuid:b650a7fd-0cb0-4396-9826-3d8731fadb09
  » http://resolver.tudelft.nl/uuid:b650a7fd-0cb0-4396-9826-3d8731fadb09
• ROBINSON RL & JACOBY RH. 1965. Better compressibility factors. Hydrocarb Process 44: 141-145.
• SAFFARI H & ZAHEDI A. 2013. A new alpha-function for the Peng-Robinson equation of state: Application to natural gas. Chinese J Chem Eng 21: 1155-1161.
• SETZMANN U & WAGNER WA.1991. New Equation of State and Tables of Thermodynamic Properties for Methane Covering the Range from the Melting Line to 625 K at Pressures up to 100 MPa. J Phys Chem Ref Data 20: 1061-1155.
• SIMON PP, HWANG CA, HOU H, HALL KR & HOLSTE JC. 1997. Burnett and pycnometric (p,Vm,T) measurements for natural gas mixtures. J Chem Thermodyn 29: 1455-1472.
• SOAVE G. 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem Eng Sci 27: 1197-1203.
• TERRON LR. 2017. Propriedades dos fluidos derivados de petróleo, 1a ed., LTC-Livros Técnicos e Científicos Editora Ltda, 524 p.
• TREBBLE MA & BISHNOI PR. 1987. Development of a new four-parameter cubic equation of state. Fluid Phase Equilibr 35: 1-18.
• TWU CH, COON JE & CUNNINGHAM JR. 1996. An approach for the application of a cubic equation of state to hydrogen-hydrocarbon systems. Ind Eng Chem Res 35: 905-910.
• U.S. EIA. 2016. Energy Information Administration. International Energy Outlook 2016 Liquid fuels [online]. Available on: https://www.eia.gov/outlooks/ieo/pdf/0484(2016).pdf Accessed 20.11.2018.
  » https://www.eia.gov/outlooks/ieo/pdf/0484(2016).pdf
• VESTFÁLOVÁ M. 2017. Possibility of derivation of energetic state variables from simplified BWR thermal state equation for water and water steam. EPJ Web Conf. 143: 1-8.
• WANG M, ZANG J, XU Q & LI K. 2011. Thermodynamic-analysis-based energy consumption minimization for natural gas liquefaction. Ind Eng Chem Res 50: 12630-12640.
• YOUNGLOVE BA & ELY JF. 1987. Thermophysical Properties of Fluids. II. Methane, Ethane, Propane, Isobutane, and Normal Butane. J Phys Chem Ref Data 16: 577-798.
• YUAN Z, CUI M, SONG R & XIE Y. 2015. Evaluation of prediction models for the physical parameters in natural gas liquefaction processes. J Nat Gas Sci Eng 27: 876-886.

Publication Dates

• Publication in this collection
  27 Nov 2020
• Date of issue
  2020

History

• Received
  24 May 2019
• Accepted
  11 July 2019