天然气水合物温度-压力-盐浓度三维相平衡曲面方程

啜世阳 王华宁,2 张宁

1.同济大学航空航天与力学学院 2.同济大学土木工程防灾国家重点实验室

赋存于土体中的固体水合物(可燃冰)对温度、压力、化学等赋存环境的改变异常敏感，升温、降压或孔隙水盐浓度增加均会导致其分解释放气体。水合物相平衡是描述水合物在平衡状态下其不同相态的温度、压力和组成等状态变量之间的关系[1]，这一关系还受多孔介质中孔径等其他因素的影响。可靠且考虑全面影响因素的水合物相平衡关系，可以正确反映水合物生成和分解条件，有助于预测水合物成藏范围、确定开采中水合物分解域等，在勘探、开采控制预测以及与水合物分解相关的工程问题中具有重要意义。本研究将尝试给出正确反映温度、压力和盐浓度这3种主要因素的影响且适用范围更广的纯水合物相平衡三维曲面方程。

目前，预测纯水合物形成/分解条件的模型主要有图解法、经验公式法、相平衡计算法以及严格的分子热力学理论模型等[2]。最后一种建立在van der Waals-Platteeuw天然气水合物统计热力学模型的基础上，其多相系统相平衡条件的理论依据涉及到水合物各相中每一组分化学势的计算，具有严格的理论推导和明确的物理意义。该理论的发展历经半个多世纪，最初由Barrer和Stuart[3]提出，van der Waals和Platteeuw[4]给出了更精确的表达，相关学者进行了完善[5-8]。

试验研究随着理论的发展逐步展开，de Roo等[9]开展了纯水中压力在11 MPa以下天然气水合物相平衡试验；Dickens和Quinby-Hunt[10]进行了海水中的相平衡试验；Maekawa等[11]、Jager等[12]及Hu等[13]将de Roo的研究工作拓展到了更高的压力水平。以上研究均将试验数据与统计热力学模型的预测进行了对比分析，并对数据进行一定的拟合处理。

此外，Clennell 等[14]和 Henry 等[15]建立了考虑单一孔径毛细效应影响的多孔介质中水合物的相平衡理论模型。

在目前的研究中，依据试验数据的相平衡拟合公式大多只考虑了在特定盐浓度下的温度-压力关系，少数研究得到了温度-压力-盐浓度三维相平衡曲面，拟合结果在精度、物理意义、使用范围等方面存在各种问题；而能够大范围较为精确地预测相平衡的统计热力学模型，又因其复杂性难以在其他工程的理论推导中直接运用。本文基于统计热力学模型和Sloan等[16]编写的CSMHYD软件[17]，通过数据拟合给出了一个描述精确、物理意义合理、适用范围广的天然气水合物温度-压力-盐浓度三维相平衡曲面方程，以方便统计热力学模型理论成果在天然气水合物开发中的应用。

1 天然气水合物三维相平衡曲面方程拟合

1.1 拟合源点与拟合函数的选取

海洋中水合物所处环境溶液(即海水)、当前普遍采用的水合物热力学抑制剂中的主要成分之一均为NaCl[18]；并且对于NaCl这种溶质的相平衡研究也最为透彻。为使本研究更有针对性且便于对研究结果的验证，针对NaCl溶液环境进一步开展了相关的相平衡研究。

根据软件CSMHYD可以找到相平衡曲面上的点作为方程的拟合点。为使点的分布更均匀、适用性更强，应尽可能覆盖较大温度、压力和溶液盐浓度范围。目前水合物热力学抑制剂中NaCl的质量分数一般在20%以内[19]，考虑到海水的稀释，工程中基本不会出现溶液环境浓度超出20%(w)的情况。同时其溶解度在25%左右，且受温度和压力影响较小，所以最终确定相平衡面中溶液NaCl质量分数X的范围是0%～25%。对于压力范围，目前已知的天然气水合物分布在水深500～4000 m，对应的压力范围是5～40 MPa[19]。但高压范围内(80～200 MPa)的数据对于天然气水合物在油气安全开发和海水淡化中也有重要意义[13]。因此，本研究相平衡面拟合中压力p范围为3～200 MPa。由此确定的温度T范围为257.8～316.8 K。据文献中拟合曲面的习惯，相平衡曲面(方程)的x、y、z轴分别为NaCl质量分数(%)、压力(MPa)、温度(K)。根据选定的范围确定拟合点的NaCl质量分数和温度，并由CSMHYD软件求出对应的压力p，共确定194个点作为拟合源点，见图1。源点在拟合范围内均有分布，且在重要位置密度更大，以保证拟合结果在此区域的准确性。

描述准确、合理且有物理含义的拟合函数形式对三维相平衡方程的拟合精度和适用性非常重要。由Clausius-Clapeyron方程可以看出，相平衡关系中临界温度T和压力p之间存在关系见式(1)。

(1)

式中：ΔHd为天然气水合物分解焓；m为压缩因子；R为气体常数。

积分可以看出临界温度与压力的对数相关，见式(2)。

(2)

式中：a、b为常数。

考虑到离子浓度对压力和温度的影响，拟合函数中认为温度与压力的对数正相关，但具体形式与上式不同。对于NaCl质量分数X和温度T之间的关系，de Roo[9]根据试验结果得到类似于温度和压力之间的关系见式(3)。

(3)

式中：c、d为常数；X′为摩尔分数，与质量分数正相关。

因此，拟合函数中认为温度与NaCl质量分数X的对数负相关。拟合源点对应的临界温度沿盐浓度方向单调递减，而自然对数函数是单调增函数，因此，拟合前将所有的点关于yOz平面对称，再向x轴正方向平移26个单位，得到一组新的拟合源点，从而方便使用自然对数函数进行拟合。拟合完毕后再进行逆向移动和对称，即可得到需要的显式函数表达式。根据上述分析，三维相平衡拟合方程形式采用分离变量形式(针对对称和平移后的拟合源点)，如式(4)所示。

T(X,p)=A(X)lnp+B(X)[lnp]3+

C(X)[lnp]5+T0D(X)

(4)

式中：A(X)=a1+a2lnX+a3[lnX]3

B(X)=b1+b2lnX+b3[lnX]3

C(X)=c1+c2lnX+c3[lnX]3

D(X)=d1+d2lnX+d3[lnX]3

ai,bi,ci,di(i=1,2,3)为常数；T0均为待拟合系数。

为提高拟合精度，假设函数中温度T与lnX和lnp的1～3次幂或者1～5次幂相关。为保证拟合函数在全域中的单调性得到满足，只保留了奇数次幂。试算发现，当幂次更高时，其对应的系数已经接近于零。拟合函数覆盖从零次到高次及其交叉项的所有可能项，可以使拟合方程对应的曲面与拟合源点更加接近。

1.2 拟合结果

T(X,p)=247.268+3.7167ln(26-X)+

0.0793175[ln(26-X)]3+

lnp{9.57275-0.547938ln(26-X)+

0.104963[ln(26-X)]3}+

[lnp]3{-0.118983+

0.0123896ln(26-X)-

0.00260223[ln(26-X)]3}+

为了进一步探求图式理论策略与阅读能力的关系，该研究采用逐步进入法（stepwise）对图式理论策略与阅读能力进行多元回归分析（见表3）。

[lnp]5{0.00388941-

0.0000814427ln(26-X)+

0.0000307206[ln(26-X)]3}

(5)

拟合优度已非常接近1。平移的单位可以由26改为其他正数，该数字决定了方程适用的NaCl的质量分数范围(本文是0%～26%，但在临近26%时与数据源点的误差较大，因此，建议在质量分数为25%以下时使用)。

图2(a)给出由拟合方程确定的曲面及其拟合源点，图2(b)为3%、9%、17%、25% 4种质量分数下的截线与拟合源点的对比。由图2可以看出，曲面/线与拟合点在整个拟合范围差别都很小。拟合方程反映了工程应用范围内天然气水合物临界相平衡处温度、压力、化学浓度的基本规律：当固定化学浓度时，临界温度随压力的增加按照自然对数的幂函数形式增加；当固定压力时，临界温度随质量分数的增加按照自然对数的幂函数形式减小。

2 相平衡方程的验证和比较

将本研究提出的相平衡方程与多篇文献中覆盖面较广的二维和三维试验数据进行对比，从而验证方程在描述实际水合物分解中的适用性。

2.1 与固定盐浓度下试验数据及拟合曲线的比较

Dickens和Quinby-Hunt[10]及Liu[20]通过试验测定了纯水和海水环境中天然气水合物的相平衡点，并给出了海水环境下的拟合曲线，见式(6)(单位已换算成与本研究一致，下同)。

(6)

此后，Brown 等[21]、Peltzer和Brewer[22]根据Dickens 海水中相平衡试验数据分别得到了另外的经验公式。

图3给出本研究结果与纯水(X=0)和海水环境(X=3.35%)试验数据以及上述3个拟合方程的对比。可见，本研究方程能够很好地满足试验结果，与 Dickens等[10]及其他文献[21-22]的拟合曲线也十分接近。但除纯水和海水的化学环境，本研究结果还可适用于0%～25%的任何盐浓度，适用范围更广。

2.2 与三维相平衡试验点的比较

de Roo等[9]在温度为261.85～285.98 K、最大压力为11.0 MPa的条件下，通过试验研究了NaCl不饱和溶液中，温度和压力对甲烷水合物形成的影响。Maekawa等[11]得到了在18 MPa以下压力范围内水合物相平衡的试验数据。这两组中低压力下实验得出的结果和本文的相平衡曲面之间的对比如图4所示。由图4可以看出，本文曲面和实验点在中低压条件时高度吻合。

Jager和Sloan[12]测量了在纯水和4种不同浓度NaCl溶液中，压力不超过70 MPa的情况下天然气水合物分解的相平衡数据。Hu等[13]的研究揭示了在高盐度和高压的极端条件下，天然气水合物的不寻常的形成过程，压力提高至约200 MPa。将这两组高压和高盐度条件下测得的数据与本文曲面对比，如图5所示，在压力较高的区域(大于50 MPa)，预测值开始出现了稍大的偏差，温度预测值均略高于试验点(15%

总体来说，本研究相平衡曲面方程的拟合函数选取更为合理，高盐度、高压区域的预测值比现有曲面精确[9,11]。

3 三维相平衡曲面部分特点的验证

3.1 化学因素对T-p曲线的影响

基于相平衡方程，图6给出不同盐浓度下的T-p曲线。由图6可见，压力较大时，各曲线之间基本满足平移关系；但在压力较小时很难直观判断。因此，下面对相平衡方程两边对X求偏导，给出温度沿盐浓度方向的变化梯度。若该偏导数在压力p方向变化很小，则意味着盐浓度改变时，其温度差值在不同压力下基本相同，表示平移关系近似成立。求偏导后的结果见式(7)。

T(X,p)x={3.7167-0.547938lnp+

0.0123896[lnp]3-

0.0000814427[lnp]5+{0.237952+

0.314889lnp-0.00780668[lnp]3+

0.0000921617[lnp]5}ln(26-X)2}/

(-26+X)

(7)

图7为式(7)对应曲面图形。可见当盐浓度较低、压力高于10 MPa时，其偏导数沿压力方向几乎不发生变化，可以认为T-p曲线满足平移关系；但是在压力较低时，其偏导沿压力方向变化较为剧烈，T-p曲线在这一区域并不能视为平移。

3.2 Clausius-Clapeyron线性行为的验证

在不考虑盐浓度影响时，天然气水合物相平衡关系满足Clausius-Clapeyron方程。Dickens等[10]、孙始财等[23]以及其他研究人员的试验都证明了这一点。为验证本文方程是否具有这一特点，绘制研究得到的1/T-X-lnp关系的曲面见图8，曲面整体接近一个平面。在不同的盐浓度下，1/T和lnp之间基本上保持线性关系，局部线性程度很好，整体有一定非线性，与前人的研究结果一致。

4 结论与展望

(1) 基于天然气水合物相平衡的统计热力学理论和CSMHYD软件，拟合得到温度-压力-NaCl盐浓度三维相平衡曲面方程。方程形式简洁、物理关系明确，适用的压力和NaCl溶液浓度区间几乎涵盖了目前工程中可能遇到的所有情况，且在描述已有的试验数据时具有很高的吻合度，特别在高压、高浓度区域比现有的拟合曲线或曲面更优，适用范围更广。

(2) 基于方程揭示了在压力高于10 MPa时，化学因素将近似导致T-p曲线平移，但压力较小时非平移所能描述；Clausius-Clapeyron线性行为在整个曲面上基本成立。

(3) 后续研究将基于更加准确的理论模型和软件(如CSMGem)以类似的方法拟合出更精确的曲面，并进一步考虑对水合物赋存于土体中时，更复杂的孔隙环境下相平衡曲面方程的修正。