基于CPA状态方程的高压含硫天然气压缩因子计算方法研究

罗召钱 李旭成 杜诚 余相君

中国石油西南油气田公司川西北气矿

天然气的压缩因子是气藏工程、采气工艺和地面集输工艺中必需的重要基础参数[1]。随着近年来我国在深层油气勘探开发领域不断取得突破，涌现出了大量的高压含硫气田。例如川西北部地区的九龙山区块茅口组气藏原始地层压力高达127 MPa，H2S质量浓度达到15 g/m3。同样，在川西北部地区的双鱼石区块等同样发现了原始地层压力高达90 MPa以上的含硫气藏[2]。与普通天然气相比，含硫天然气中的H2S含量将对天然气压缩因子计算产生重要影响。

目前，国内外学者对压缩因子计算主要是依据经验公式和状态方程开展[3]。其中，常用的经验公式包括Dranchuk-Purvis-Pobinson方法、Brill-Beggs方法、Pedlich-Kwong方法；这些经验公式虽然在一定的压力、温度范围内具有较高的精度，但难以扩展应用到比热容、焓、熵等其他热物性参数和相平衡的计算[4]。而状态方程不仅可以对压缩因子进行计算，还可以对天然气的其他物性参数和相平衡进行准确的计算。

国内外一般采用SRK、PR等立方型状态方程计算天然气的压缩因子。但是，SRK、PR等立方型状态方程主要是针对非极性体系推导的，而H2S是强极性分子，特别是高压条件下分子之间的间距更小，极性分子之间的作用力更强，传统状态方程计算压缩因子存在一定误差。CPA(Cubic-Plus-Association)状态方程在传统立方型状态方程的基础上，引入考虑极性分子间相互作用的缔合项，克服了传统立方型状态方程难以描述含极性物质物性参数的不足[5]。因此，基于CPA状态方程有望对高压含硫天然气的物性参数进行更加准确的计算。

针对上述问题，本文基于文献数据对PR、SRK、CPA 3种状态方程计算高压含硫天然气压缩因子的精度进行了评价，修正了CPA状态方程中的H2S与CH4的二元交互作用系数。在此基础上，采用法国ST抗硫高压PVT仪验证了压力为70～131.1 MPa范围内基于CPA状态方程计算含硫天然气压缩因子的精度。

1 压缩因子计算状态方程

1.1 CPA状态方程

目前，用于计算天然气压缩因子的状态方程多达上百种，PR、SRK状态方程是在油气行业应用较为广泛的状态方程。然而，SRK、PR等立方型状态方程只考虑了非极性分子之间的色散力、诱导力，没有考虑H2S等极性分子之间的缔合作用。在压力低于20 MPa时，分子之间的间距较大，缔合作用较弱，精度较高。然而，当压力高达70 MPa以上时，H2S分子之间的间距减小，分子之间的缔合作用增强，导致基于传统的立方型状态方程并不适用于含H2S天然气的物性参数计算。

丹麦科技大学Kontogeorgis教授在传统立方型状态方程的基础上，通过引入考虑极性分子间缔合作用的缔合项，提出了CPA状态方程。CPA状态方程的基本形式如式(1)所示[5-6]。

(1)

式中：p为压力，kPa；T为温度，K；R为气体常数；a为能量参数，(kPa·m6)/mol2，采用式(2)计算；b为体积参数，m3/mol；v为摩尔体积，m3/mol；g为分子的半径分布函数，采用式(3)计算，无量纲；xi为混合物中组分i的摩尔分数，无量纲；Ai代表分子i上的活性缔合点位A；XAi为组分i中未缔合的活性点位A的摩尔分数，采用式(4)计算。式(1)中等号右边的3项分别为斥力项、引力项和缔合项。

(2)

式中：Tr=T/Tc，Tc为临界温度，K；κ为修正系数；a0为能量参数。

(3)

式中：η=b/4v。

(4)

式中：NCA为极性组分的数量；ρ为混合物摩尔密度，mol/m3;ΔAiBj为分子i上的缔合点位A与分子j上缔合点位B之间的缔合强度，采用式(5)计算：

(5)

(6)

式中：εAiBj为活性缔合点位Ai与Bj之间的缔合能量，(kPa·m3)/mol；βAiBj为活性缔合点位Ai与Bj之间的缔合体积，无量纲。

对于混合物，a和b采用van der Waals混合规则计算；εAiBj和βAiBj采用CR-1混合规则计算，如式(7)～式(10)所示：

(7)

(8)

(9)

(10)

式中：kij为二元交互作用系数。

CPA状态方程采用a0、b、κ、εAiBj、βAiBj共5个参数对每一种组分进行表征。对于非极性组分(甲烷、乙烷等)，εAiBj=0、βAiBj=0，所以CPA状态方程退化为SRK状态方程。因此，在计算高压含硫天然气的压缩因子时，关键在于获取H2S分子的缔合参数。

1.2 H2S分子的缔合参数

H2S通常被认为是非缔合化合物，但Cabaleiro-Lago 等[7]研究了萘-H2S的相互作用，认为H2S存在较弱的自缔合倾向。Ioannis Tsivintzelis等[8]对比了H2S与H2O、CO2之间的缔合模型，发现对于与水和碳氢化合物的混合物，采用3B缔合方案将H2S视为自缔合分子并假设CO2仅具有1个质子受体位点的溶剂化时，CPA状态方程的计算值更接近实验值。因此，对于H2S采用3B缔合方案，即1个电子接受点位、2个电子贡献点位，缔合参数β=0.232 9，ε=3.781(Pa·m3)/mol[9]。

2 高压含硫天然气压缩因子计算精度评价

综合考虑压力和H2S含量两个因素，根据SY/T 6581-2012 《高压油气井测试工艺技术规程》和GB/T 26979-2011《天然气藏分类》将气井分为7类，分别对SRK、PR和CPA方程进行适应性评价，如表1所示。

基于文献中公开发表的数据，调研了压力为3.72～97.58 MPa、φ(H2S)为0%～70.03%的天然气压缩因子数据共154组，包括郭绪强等[10]测定的压力为20～100 MPa的压缩因子数据；普光气田普光6井压力为8.5～50 MPa、φ(H2S)为14.99%的数据[1]；河坝1井压力为15～43 MPa、φ(H2S)为0.675%的数据[11]；罗家寨气藏飞仙关组鲕滩气藏压力为10～40 MPa、φ(H2S)为8.36%的数据[12]；梁光川等[13]测定的压力为7～49 MPa的天然气压缩因子；Elsharkawy A M等[14]实验测试的φ(H2S)为6.8%～70.03%、压力为6.89～37.13 MPa的压缩因子；Zhou等[15]测试的压力为3.72～18.88 MPa的天然气压缩因子。将上述实验数据按照含硫量和气藏压力分类，结果如表2所示。

表1 按压力、硫含量分类的气井Table 1 Gas wells classified by pressure and sulfur contents类别φ(H2S)/%p/MPa中低压低含硫气藏<0.3中低压中高含硫气藏0.3≤φ<10中低压特高含硫气藏≥10<35高压低含硫气藏<0.3高压中高含硫气藏0.3≤φ<10高压特高含硫气藏≥1035≤p<70超高压气藏0≤φ≤100≥70

表2 按压力、硫含量分类的压缩因子实验数据Table 2 Experimental compressibility factor data classified by pressure and sulfur contents类别压缩因子实验数据数量/个中低压低含硫气藏15中低压中高含硫气藏14中低压特高含硫气藏42高压低含硫气藏2高压中高含硫气藏7高压特高含硫气藏14超高压气藏60

CPA状态方程中并未提供H2S与CH4、CO2等常见组分之间的二元交互作用系数，这会影响到压缩因子的计算精度。为此，基于文献中发表的压缩因子，

表3 H2S与CH4、CO2分子间的二元交互作用系数[16]Table 3 Binary interaction parameters between H2S and CH4, CO2 molecules极性组分CH4H2SCO2H2O-0.039 80.081 6-0.055 9CH40.000 00.088 80.095 6H2S0.088 80.000 00.115 0CO20.095 60.115 00.000 0

修正了CPA状态方程中H2S与CH4、CO2间的交互作用系数，结果如表3所示。根据表3将上述实验数据与SRK、PR和CPA状态方程的计算结果进行了对比，结果如图1～图7所示。

定义压缩因子的文献数据与状态方程计算值的相对偏差θ如式(11)所示，SRK、PR、CPA状态方程的压缩因子与文献值之间的偏差如表4和表5所示。

(11)

表4 状态方程压缩因子计算值与文献值的平均相对偏差Table 4 Average relative deviations between calculated compressibility factors and literature data气藏分类平均相对偏差/%SRKPRCPA中低压低含硫气藏2.92-2.572.86中低压中高含硫气藏4.15-1.143.05中低压特高含硫气藏4.84-0.602.23高压低含硫气藏1.65-6.400.24高压中高含硫气藏3.78-2.692.39高压特高含硫气藏6.32-1.541.46超高压气藏-2.70-10.68-2.63

表5 状态方程压缩因子计算值与文献值的最大相对偏差Table 5 Maximum relative deviations between calculated compressibility factors and literature data 气藏分类最大相对偏差/%SRKPRCPA中低压低含硫气藏4.39-4.552.86中低压中高含硫气藏5.13-4.224.96中低压特高含硫气藏7.20-4.667.16高压低含硫气藏2.03-6.890.48高压中高含硫气藏4.68-4.233.60高压特高含硫气藏7.63-2.522.36超高压气藏-3.39-11.213.07

综合分析表4和表5可得，对于中低压含硫天然气(p<35 MPa)，采用PR方程计算精度最高，计算得到平均相对偏差为1.12%；对于高压及超高压含硫天然气(p>35 MPa)，采用CPA方程计算精度最高，计算得到平均相对偏差为-1.46%。这是因为在压力低于35 MPa时，H2S分子之间的缔合作用较弱，缔合作用对于压缩因子的影响可以忽略，而PR状态方程是在SRK状态方程的基础上改进而来的。因此，PR状态方程计算的压缩因子比基于SRK状态方程发展而来的CPA状态方程的精度更高。然而，随着压力的升高，H2S分子之间的缔合作用增强，H2S分子对于压缩因子的影响逐渐凸显，因此在压力高于35 MPa时，CPA状态方程的精度更高。

3 超高压含硫天然气压缩因子实验测试与计算

3.1 超高压含硫天然气压缩因子的实验测试

为了进一步验证采用CPA状态方程计算超高压含硫天然气压缩因子的精度，采用法国ST抗硫高压PVT仪测试了川西北部地区4组超高压含硫天然气在压力为70.00～131.10 MPa内的压缩因子。ST抗硫高压PVT仪流程如图8所示，主要包括高压计量泵、储样容器、PVT筒和恒温箱；最大工作压力为150 MPa，最大工作温度200 ℃，压力控制精度可到0.01 MPa，温度控制精度到0.1 ℃，可完成φ(H2S)≤20%、φ(CO2)≤50%的天然气的PVT测试。

表6 超高压高温含硫天然气压缩因子实验数据Table 6 Experimental data of ultra-high pressure and high temperature sulfur-containing natural gas compression factor气样压力范围/MPaφ(H2S)/%数据数量/个L004-X1井75.10～131.100.878 339L016-H1井72.00～105.000.000 133ST1井73.00～122.880.025 033SY001-1井70.00～100.000.356 433

实验共测得压力为70.00～131.10 MPa，温度为100～157 ℃，φ(H2S)为0.000 1%～0.878 3%共138组超高压高温含硫天然气压缩因子实验数据，见表6。

3.2 状态方程压缩因子计算值与实验值对比分析

将实验值与SRK、PR及CPA状态方程的计算值进行对比，如图9所示。

由图9分析可见，对于川西北地区超高压含硫天然气压缩因子的计算，CPA状态方程精度最高，平均相对偏差为-1.50%，最大相对偏差为-2.99%；SRK状态方程其次，平均相对偏差为-1.72%，最大相对偏差为-3.65%；PR的精度最低，平均相对偏差为-9.18%，最大相对偏差为-11.36%。因此，本研究修正的CPA状态方程能够准确地计算超高压含硫天然气的压缩因子，该方法还可以进一步推广应用于超高压含硫天然气的气液相平衡计算，为气藏工程设计、地面集输工艺设计中相关物性参数的计算提供依据。

4 结论

(1) 对于中低压含硫天然气(p<35 MPa)，采用PR方程精度更高，平均相对偏差为1.12%；对于高压及超高压含硫天然气(p>35 MPa)，采用CPA方程精度更高，平均相对偏差为-1.46%。

(2) 开展了川西北地区超压含硫天然气压缩因子实验测试，将测试的138组实验数据与SRK、PR、CPA状态方程计算值进行了对比。结果表明，对于超高压含硫天然气，CPA方程具有最高精度，且平均相对偏差为-1.50%，最大相对偏差为-2.99%。

(3) 在高压下，气体分子之间的间距减小，分子间作用力的影响不可忽视，而CPA方程考虑了H2S等极性分子之间的相互缔合关系，因此，对于天然气参数计算具有更高的精度。