石油与天然气化工  2016, Vol. 45 Issue (4): 38-41, 50
天然气体系三相平衡计算交互作用因子Kij的取值
李浩勇 1, 王登海 2, 耿树博 1, 郑欣 2, 赵丹 1, 王文珍 1     
1. 西安石油大学;
2. 长庆科技工程有限公司
收稿日期：2016-01-25
作者简介：李浩勇(1992-)，男，硕士研究生。现就读于西安石油大学，研究方向为天然气体系热力学计算及LNG热力学条件与天然气体系脱硫脱碳的研究。E-mail：lihaoyong10@163.com.
摘要：为解决LNG产业链中的重烃析出问题，研究建立了天然气体系三相相平衡理论模型并自主研发了热力学计算软件。计算模型中的Kij代表纯物质变成混合物时所有的物理性质变化，是温度、组成以及组分性质的函数。目前对于物质混合过程的物理化学变化没有全面的认识，所以Kij没有具体的表达式。计算了天然气混合体系中不同Kij的相平衡计算，以选取合适的Kij值来提高模型计算结果的准确性。
关键词：液化天然气    交互作用因子    相平衡    逸度计算    
Calculation of interaction parameter on phase equilibrium in natural gas
Li Haoyong1 , Wang Denghai2 , Gen Shubo1 , Zheng Xin2 , Zhao Dan1 , Wang Wenzhen1     
1. Xi'an Shiyou University, Xi'an 710065, China;
2. Changqing Science and Technology Engineering Co., Ltd., Xi'an 710065, China
Abstract: To solve the heavy hydrocarbon precipitation problems in LNG industry, the theoretical model of three phase equilibrium is established and the thermodynamic calculation software is developed independently. The interaction parameter Kij in calculation model represents all the physical properties changing when the pure substance becomes mixture, and it's the function of T, x and component properties. It can't be calculated by a formula as physical chemistry changes of the mixing process have not been understood comprehensively now. The phase equilibrium in different Kij is calculated in order to select the appropriate Kij to improve the accuracy of the model calculation results.
Key Words: liquefied natural gas    interaction parameter    phase equilibrium    fugacity calculation    

我国天然气消费量大，东西部天然气分布与消费极不平衡，液化天然气技术成为我国中短距离运输天然气的主要方法[1-3]。在液化流程中，天然气在许多设备中处于低温状态，其重烃等组分在低温下会形成固相，堵塞管道、阀门，给生产造成事故和危害。所以，有效地脱除天然气中的重烃是天然气液化流程中非常重要的技术环节[4-11]。本文通过自主研发的热力学计算软件，研究建立了天然气体系三相相平衡理论模型并进行了计算，为天然气的液化工艺设计提供了可靠的数据支持和基础理论依据。

建立完善的相平衡计算体系，首先要保证计算结果的准确性和精度。在本模型中有很多影响计算精度的因素，如状态方程、混合规则以及交互作用因子[12-15]。状态方程和混合规则都可以根据计算目标体系的性质来选取，以便将误差减至最小，而交互作用因子Kij代表了多种纯物质变成混合体系时所有的物理性质的变化。由于目前对混合过程的物理化学机理认识还不够充分，没有一个普遍通用的Kij函数具体表达形式[16-17]。本文根据计算来确定最适合的Kij参数值，通过计算不同Kij的液相逸度并与Aspen计算结果进行分析比较，最终确定本模型Kij的最优值。

1 本模型计算的原理

本模型气液两相平衡用PR方程计算：

$ p=\frac{RT}{V-b}-\frac{a\left( T \right)}{V\left( V+b \right)+b(V-b)~} $ (1)

式中：$ a=\sum\limits_{i}{\sum\limits_{j}{{{x}_{i}}{{x}_{j}}{{a}_{ij}}}} $

$ \begin{align} &b=\left( 1-{{K}_{ij}} \right){{\left( {{\alpha }_{i}}{{\alpha }_{j}}{{\beta }_{i}}{{\beta }_{j}} \right)}^{0.5}} \\ &{{\alpha }_{i}}=0.457\ 24{{R}^{2}}T_{ci}^{2}/{{p}_{ci}} \\ &{{\beta }_{i}}=[1+(0.374\ 64+1.542\ 26{{\omega }_{i}}-\\ &\ \ \ \ \ \ \ \ \ \ 0.269\ 92\omega _{i}^{2})\left( 1-T_{ri}^{0.5} \right){{]}^{2}} \\ \end{align} $
$ \begin{align} &b=\sum\limits_{i}{{{x}_{i}}{{b}_{i}}} \\ &~{{b}_{I}}=0.077\ 8R{{T}_{ci}}/{{p}_{ci}} \\ \end{align} $

Kij为交互作用因子，将在第2节讨论。

令A=ap/(RT)2, B=bp/RT, Z=pV/RT，则组分i的液相逸度系数φil用式(2)计算：

$ \begin{align} &\text{ln}\varphi _{i}^{l}=\frac{{{B}_{i}}}{B}(Z-1)-\text{ln}\left( Z-B \right)+\ \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{A}{2.828}\left[\frac{{{B}_{i}}}{B}-\frac{2}{a}\sum{{{x}_{ij}}{{a}_{ij}}} \right]\text{ln}\left[\frac{Z+2.414B}{Z-2.414B} \right] \\ \end{align} $ (2)

所以组分i的液相逸度fil=xiφilp，其中，xi为组分i的摩尔分数，p为体系压力。

本模型将固相的析出组分作为纯固体物质处理，计算出固体标准逸度fios，进而计算天然气中该组分的固相逸度fis，对于固相其fis≈fios，其中标准逸度由以下热力学关系式(式3)推出。

$ f_{io}^{s}=f_{io}^{l}\text{exp}\left[-\frac{\Delta H_{i}^{f}}{RT}\left( 1-\frac{T}{T_{i}^{f}} \right)+\ \frac{1}{RT}\int _{{{T}^{i}}}^{{{T}^{f}}}\left( C_{pi}^{l}-C_{pi}^{s} \right)dT-\ \frac{1}{R}\int _{{{T}^{i}}}^{{{T}^{f}}}\frac{\left( C_{pi}^{l}-C_{pi}^{s} \right)}{T}dT \right] $ (3)

式中：fiol为组分i的液相标准态逸度，fil=xi fiol；(Cpil-Cpis)为恒压热容差，ΔCpi=Cpil-Cpis=b1Mi+b2MiT，b1和b2为常数，取近似值为b1=0.303 3，b2=-4.635×10-4；ΔHif为组分i的熔化焓；Tif为组分i的熔点。

2 交互作用因子Kij的取值
2.1 不同条件下Kij的取值

文献资料显示[17-18]，多组分混合体系中Kij的取值一般在0到0.1之间。本文在0到0.1之间对Kij取值，计算了混合体系的液相逸度系数。计算体系组成及性质如表 1所示[19]。

表 1    计算体系组成及性质 Table 1    Composition and property of calculation system

由计算数据分析得出，Kij的取值和烃类中碳原子的数量有关系。以碳四和碳五为分界线，碳四及以下轻烃和碳五及以上重烃的Kij取值规律不同。

2.1.1 碳四及碳四以下轻烃的Kij取值

碳四及碳四以下轻烃在Kij=0时的计算结果非常接近Aspen数据，误差很小：当Kij值越大时，误差越大。图 1~图 5为碳四以下的轻烃在Kij=0时的液相逸度系数φ和Aspen数值对比。由图 1~图 5中可以看出，Kij=0时，除了甲烷的误差在温度大于200 K后略微明显，其余组分液相逸度系数计算值和Aspen计算值均比较接近。

图 1     甲烷的液相逸度系数计算 Figure 1     Calculation of methane liquid phase fugacity coefficient

图 2     乙烷的液相逸度系数计算 Figure 2     Calculation of ethane liquid phase fugacity coefficient

图 3     丙烷的液相逸度系数计算 Figure 3     Calculation of propane liquid phase fugacity coefficient

图 4     异丁烷的液相逸度系数计算 Figure 4     Calculation of isobutane liquid phase fugacity coefficient

图 5     正丁烷的液相逸度系数计算 Figure 5     Calculation of butane liquid phase fugacity coefficient

2.1.2 碳五及碳五以上重烃的Kij取值

碳五及碳五以上重烃的Kij取值规律与温度有关。本文选取重烃里比较重要的两个组分环己烷和苯为例来说明重烃的Kij取值与温度的关系。表 2为本模型计算的两种重烃的液相逸度系数与Aspen结果对比误差较小时的Kij取值情况。

表 2    环己烷和苯的液相逸度系数 Table 2    Liquid phase fugacity coefficient of cyclohexane and benzene

根据表 2的数据绘制了重烃的Kij取值和温度的关系曲线，见图 6。重烃的Kij随温度的变化规律可以分为两部分：当温度小于160 K时，Kij随温度的变化为线性关系；当温度大于160 K时，Kij随温度的变化为高次多项式，关系式如式(4)所示。

图 6     重烃Kij值与温度的关系 Figure 6     Relation of Kij value of heavy hydrocarbon with temperature

本模型重点解决天然气中的重烃析出问题，所以对重烃的计算尤为重要。利用式(4)作为Kij的取值，计算了苯和环己烷的液相逸度系数，见图 7和图 8。由图可以看出，当Kij按式(4)计算时，重烃的相平衡计算结果和Aspen的计算结果很接近，误差很小。

$ {{K}_{ij}}=\left\{ \begin{align} &-0.000\ 125T+0.06\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ T<160\text{ K} \\ &2.51\times {{10}^{-8}}{{T}^{3}}-1.574\ 2\times {{10}^{-5}}{{T}^{2}}+0.002\ 9T-0.128\ 4\ \ \ \ \ \ \ T\ge 160\ \text{K} \\ \end{align} \right.~ $ (4)
图 7     不同温度下苯的液相逸度系数 Figure 7     Liquid phase fugacity coefficient of benzene at different temperature

图 8     不同温度下环己烷的液相逸度系数 Figure 8     Liquid phase fugacity coefficient of cyclohexane at different temperature

2.2 压力对Kij值的影响

本研究利用5.0 MPa下的Kij公式计算了3.0 MPa和7.0 MPa下的相平衡数据。表 3和表 4列出了丙烷和苯的相平衡数据，丙烷代表碳四及以下的轻烃，苯代表碳五及以上的重烃。

表 3    3.0 MPa下苯和环己烷的液相逸度系数 Table 3    Liquid phase fugacity coefficient of cyclohexane and benzene at 3.0 MPa

表 4    7.0 MPa下苯和环己烷的液相逸度系数 Table 4    Liquid phase fugacity coefficient of cyclohexane and benzene at 7.0 MPa

由表 3和表 4可以看出，两个压力下的计算结果随Kij的变化规律和5.0 MPa下的变化规律相同。代表碳四及以下轻烃的丙烷在Kij=0时计算误差小，代表碳五及以上重烃苯的Kij取值和5.0 MPa下的取值一样。在此计算模型中，压力对交互作用因子的取值几乎没有影响。

3 结论

本文利用自主研发的热力学软件对天然气体系进行了相平衡计算，为了找到适合模型的Kij值，分别取不同的Kij值计算了不同温度压力下的液相逸度系数，并与Aspen的计算结果进行对比。结果显示，对于碳四及以下的轻烃，当Kij=0时，其计算结果准确度较高；对于碳五及以上的重烃，Kij值为温度的函数，当Kij按式(4)取值时，计算结果准确度较高。

参考文献
[1]
黄志光, 汪荣顺, 石玉美, 等. 小型天然气液化装置的研制现状与前景[J]. 低温工程, 2002, 130(6): 59-62.
[2]
邢云, 刘淼儿. 中国液化天然气产业现状及前景分析[J]. 天然气工业, 2009, 29(1): 120-123.
[3]
徐文渊. 小型液化天然气生产装置[J]. 石油与天然气化工, 2005, 34(3): 161-164.
[4]
HOSSEIN BEIGI H R, EDALAT M. A two parameters equation of state for phase equilibrium behaviour reservoir and liquified natural gas[J]. Society of Petroleum Engineers, 2005(6): 1-22.
[5]
赵敏, 厉彦忠. C3/ MRC液化流程中原料气成分及制冷剂组分匹配[J]. 化工学报, 2009, 60(增刊1): 50-57.
[6]
SERDYUKOV, LENTRANSGAZ, KHODORKOV. Mini plants at pressure reduction stations in Russia[J]. LNG Journal, 2002(4): 6-9.
[7]
李廷勋, 郭开华. LNG安全规范现状[J]. 天然气工业, 2008, 27(6): 1-3.
[8]
皮艳慧, 廖柯熹, 孙欧阳. 天然气水合物生成物条件预测模型及适用性评价[J]. 天然气与石油, 2012, 30(6): 16-18.
[9]
陈坚, 黄志明, 包永忠. 烷烃混合物二元相互作用参数的新型关联式[J]. 化工学报, 2006, 57(7): 1508-1512.
[10]
姚军辉, 石勇. 多孔介质中天然气水合物分解热力学模型研究[J]. 石油与天然气化工, 2009, 38(1): 11-14.
[11]
曾凡平, 田广新. 液化天然气工厂重烃脱除工艺方案比选[J]. 石油与天然气化工, 2015, 44(1): 35-38.
[12]
陈振华, 曹堃, 姚臻. 基于PR状态方程的二元体系临界性质计算中不同混合规则对比[J]. 高校化学工程学报, 2008, 22(3): 365-370.
[13]
JAMES M P, FRENKEL D. Numerical prediction of the melting curve of n-octane[J]. Journal of Chemical Physics, 1999, 111(4): 1500-1510.
[14]
杨帆, 陈保东, 姜文全, 等. PR方程在天然气热物性计算中的应用研究[J]. 当代化工, 2013, 42(5): 649-656.
[15]
于志家, 杨筱恬. 用SRK方程与PR方程求算双组份混合气体热力学性质[J]. 化工高等教育, 2013, 30(1): 59-62.
[16]
SQUILLACOTE M, SHERIDAN R S, CHAPMAN O L. Spectroscopic detection of the twist-boat conformation of cyclohexane[J]. J Am Chem Soc, 1975, 97(11): 3244-3246. DOI:10.1021/ja00844a068
[17]
许国华, 黄永华, 张鹏. 流体混合物热力学性质计算的混合法则新进展[J]. 制冷技术, 2002, 36(8): 31-36.
[18]
李延勋, 刘志刚. 基于PR方程的气液相平衡计算[J]. 广东工业大学学报, 1999, 16(1): 29-34.
[19]
王文珍, 王登海, 郑欣. 天然气中低碳烷烃的汽-液平衡研究[J]. 广州化工, 2014, 42(16): 13-16. DOI:10.3969/j.issn.1001-9677.2014.16.005
X