A Thermodynamic Model for Predicting 3D Phase Equilibrium Surface of Natural Gas Hydrates Considering Salt Concentration
-
摘要: 准确预测含盐体系中天然气水合物的热力学稳定性,对预测水合物成藏范围、确定水合物分解域等具有重要意义。因此,综合考虑水合物晶穴非球形特征和含盐体系中电解质相互作用,建立考虑盐浓度的天然气水合物相平衡热力学模型,并与实验数据进行对比验证。研究结果表明,建立的三维相平衡曲面能有效预测天然气水合物热力学稳定性,纯水条件下的温度绝对平均相对偏差仅为0.08,且含盐条件下不超过0.15。在氯盐摩尔分数大于0.02、压力高于20 MPa时,化学因素导致p-T曲线发生平移。在压力较小时,平衡温度沿盐浓度方向的梯度剧烈变化,p-T曲线不再具有平移特性。同时,摩尔分数相同时,AlCl3比其他氯盐对CH4水合物的抑制作用更强。CH4水合物的lnp-X-1/T三维曲面局部表现出较好的Clausius-Clapeyron线性行为,曲面整体具有非线性特征,且氯盐电解质对CH4水合物的抑制效果越好,非线性特征越强。
-
关键词:
- 天然气水合物 /
- 热力学模型 /
- 相平衡 /
- Clausius-Clapeyron线性特征
Abstract: Accurate prediction of the thermodynamic stability of gas hydrates in salt-bearing systems is of great significance for predicting the accumulation range and determining the decomposition domain of hydrate. Therefore, considering the non-spherical characteristics of hydrate crystal cavity and the electrolyte interaction in the salt-containing system, a thermodynamic model of gas hydrate phase equilibrium considering salt concentration was established, and compared with the experimental data. The results show that the established three-dimensional phase equilibrium surface can effectively predict the thermodynamic stability of NGH, and the absolute mean relative deviation of temperature is only 0.08 in pure water, and no more than 0.15 in salt. When the molar fraction of chloride is greater than 0.02 and the pressure is higher than 20 MPa, the chemical factors cause the p-T curve to shift. When the pressure is small, the equilibrium temperature changes dramatically along the gradient of salt concentration, and the p-T curve no longer has translation characteristics. At the same time, AlCl3 has stronger inhibitory effect on CH4 hydrate than other chlorine salts when the molar fraction is the same. The lnp-X-1/T three-dimensional surface of CH4 hydrate exhibits better Clausius-Clapeyron linear behavior locally, and the surface as a whole has certain nonlinear characteristics. Moreover, the better the inhibition effect of chloride electrolyte on CH4 hydrate, the stronger the nonlinear characteristics. -
0. 引言
天然气水合物作为一种新型清洁能源,被视为21世纪最具潜力的接替能源之一[1-2]。天然气水合物对温度、压力、化学等赋存环境的改变异常敏感,温度、压力和盐浓度的轻微改变,都可能会导致水合物生成与分解[3-5]。准确了解含盐体系中天然气水合物的热力学稳定性有助于预测水合物成藏范围、确定开采过程中水合物分解域等,在勘探、开采控制预测以及与水合物形成与分解相关的工程问题中具有非常重要的意义[6-9]。
对含盐体系中气体水合物相平衡热力学研究理论模型主要为vdW-P模型和Chen-Guo模型[10]。Javanmardi[11]、Mohammadi[12]、Bandyopadhyay[13-14]、Haghighi[15]和Hsieh[16]等人基于vdW-P模型,分别使用不同的状态方程和活度系数模型,建立了天然气水合物相平衡热力学模型预测含有不同浓度NaCl、KCl、CaCl2、MgCl2水溶液体系中气体水合物相平衡。Ma[17]和Li等人[18]则基于Chen-Guo热力学模型,结合Pitzer-Debye-Hückel方程和N-NRTL-NRF活度系数模型,采用Valderrama-Patel-Teja 状态方程研究了含电解质水溶液体系中气体水合物的相行为。尽管vdW-P模型和Chen-Guo模型在一定温压范围内能准确预测含盐体系中天然气水合物的相平衡条件,但当体系中加入多种盐类电解质时,由于未考虑电解质之间的相互作用,导致预测精度不足。另外,vdW-P模型计算Langmuir吸附常数时假设球对称势,导致预测结果误差相对较大;而Chen-Guo模型中的水活度计算对状态方程非常敏感,状态方程的选取对预测结果影响非常大。并且,Chen-Guo模型中计算Langmuir吸附常数用到二元相互作用参数依赖于实验VLE数据,导致适用的电解质浓度和温压范围有限。最为重要的是,现有大多数研究只考虑了在特定盐浓度下的温度-压力关系,少数研究得到了温度-压力-盐浓度三维相平衡曲面,且其计算结果在精度、物理意义和使用范围等方面存在各种问题。
因此,考虑水合物晶穴的非球形特征,采用三层球模型来描述水分子与客体分子间的相互作用修正vdW-P模型。同时,考虑水溶液中电解质之间的相互作用,利用Pitzer-Debye-Hückel方程和N-NRTL-NRF模型确定水的活度;进一步使用可以在更宽的温压范围内准确地描述不同物系的p-V-T关系的BWRS状态方程计算气相逸度,建立温度-压力-盐浓度三维相平衡曲面,实现含盐体系中的天然气水合物热力学稳定性的准确预测。
1. 热力学模型
预测气体水合物生成条件的热力学模型的建立是依据水在水合物相和富水相中的化学位相等这一平衡条件[19],即
$$ \mu^{H} = \mu^{W}$$ (1) 式中,μH和μW分别表示水在水合物相和富水相的化学位。若以完全空的水合物相β(晶格空腔未被水分子占据)的化学位μβ为基准,其平衡条件表示为:
$$ \left\{\begin{array}{l} \Delta \mu^{H}=\Delta \mu^{W} \\ \Delta \mu^{H}=\mu^{\beta}-\mu^{H} \\ \Delta \mu^{W}=\mu^{\beta}-\mu^{W} \end{array}\right. $$ (2) 1.1 水合物模型
1.1.1 水在水合物相中的化学电位差
假设水合物晶腔为球形,且每个晶腔最多只能容纳一个气体分子,不同晶腔的客体分子间没有相互作用,水分子对水合物自由能的贡献与其所包容的客体分子大小及种类无关,并且电解质均不存在于水合物相。则水在水合物空腔内被客体分子占据前后的化学位偏差ΔμH为:
$$ \Delta\mu^{\mathrm{\mathit{H}}}=-RT\sum_i^2v_i\ln\left(1-\theta_i\right) $$ (3) 式中,
$ {\upsilon }_{i} $ 为每个水分子i型晶腔的数量;$ {\theta }_{i} $ 表示i型晶腔被客体分子所占据的比率,如下:$$ \theta_{i}=\frac{C_{i} f_{g}(T, p)}{1+{\sum\limits_i} C_{i} f_{g}(T, p)} $$ (4) 式中,
$ {C}_{i} $ 为客体分子在i类晶腔中的气体吸附常数;$ {f}_{g}(T,p) $ 为温度T、压力p的气体逸度。Van der Waal和Platteeuw通过假设球对称势来计算,如下:
$$ C_{i}(T)=\frac{4 \text{π}}{k T} \int_{0}^{R_{c}} \exp \left(-\frac{\omega(r)}{k T}\right) r^{2} d r $$ (5) 式中,k表示Boltzman常数;
$ R $ 表示晶腔半径;$ \omega \left(r\right) $ 表示水合物晶格空腔中客体分子与构成晶腔的水分子间的势能之和。实际上晶腔并非球形,即水分子与晶腔中心距离不等,采用三层球模型来描述水分子与客体分子间相互作用,其晶腔总势能
$ \omega \left(r\right) $ 由各层壳势能$ {\omega }_{i}\left(r\right) $ 加和得到[20],如下所示:$$ \omega(r)=\omega_{1}(r)+\omega_{2}(r)+\omega_{3}(r) $$ (6) 利用Kihara的势能模型计算晶腔总势能[21],表达式如下:
$$ \omega(r)=2 {\textit{z}} \epsilon \left[\frac{\sigma^{12}}{R_{c}^{11} r}\left(\delta^{10}+\frac{a}{R_{c}} \delta^{11}\right)-\frac{\sigma^{6}}{R_{c}^{5} r}\left(\delta^{4}+\frac{a}{R_{c}} \delta^{5}\right)\right] $$ (7) 式中,z为配位数;a为分子核半径,Å;
$ \epsilon $ 为能量参数,J;$ \delta $ 为势能为零时的分子核间距,Å。部分气体的Kihara势能参数和水合物分子几何参数如表1、表2所示。$ \delta $ N的表达式为:表 1 气体的Kihara势能参数气体类型 a /Å $ \sigma $ /Å $ \epsilon /k $ /Å CH4 0.3834 3.1650 154.54 C2H6 0.6760 3.1383 190.80 C3H8 0.8340 3.1440 194.55 N2 0.5290 0.2569 150.03 H2S 0.4920 3.1774 198.53 CO2 0.1773 2.9605 170.97 表 2 水合物分子几何参数[22]水合物
结构腔型 分子核
间距/Å腔体/
单元壳配位数 H2O分子/
单元壳I 512 3.95 2 20 46 51262 4.33 6 24 II 512 3.91 16 20 136 51264 4.73 8 28 $$ \delta^{N}=\frac{2}{N}\left[\left(1-\frac{r}{R_{c}}-\frac{a}{R_{c}}\right)^{-N}-\left(1+\frac{r}{R_{c}}-\frac{a}{R_{{c}}}\right)^{-N}\right] $$ (8) 式中,N表示指数,分别取值4、5、10或11。
考虑晶腔非球形,引入扰动因子
$ {Q}^{*} $ 来校正球形晶腔的Langmuir吸附常数$ {C}_{i}^{*} $ ,即$$ C_{i}=Q^* C_{i}^{*} $$ (9) 球形晶腔的Langmuir吸附常数
$ {C}_{i}^{*} $ 和校正因子$ {Q}^{*} $ 的计算公式如下:$$ \left\{\begin{array}{l}Q^*=\mathrm{exp}\left[-a_0\left(\omega\cdot\dfrac{\sigma}{R_c-a}\cdot\dfrac{\varepsilon}{kT_0}\right)^{n_0}\right] \\ C_i^*=\dfrac{4\text{π}}{kT}\int_0^{R_c}\mathrm{exp}\left[-\dfrac{\omega_1\left(r\right)+\omega_2\left(r\right)+\omega_3\left(r\right)}{kT}\right]r^2\mathrm{d}r\end{array}\right. $$ (10) 式中,
$ {n}_{0} $ 、$ {a}_{0} $ 表示晶腔的特征常数,$ \omega $ 表示客体分子的偏心因子,$ {T}_{0} $ 基准温度,一般取值为$ {T}_{0}=273.15 \,\, {\mathrm{K}} $ 。公式(10)中的相关参数见表3。表 3 三层球模型相关参数结构及腔体类型 第一层 第二层 第三层 $ {n}_{0} $ $ {a}_{0} $ $ {R}_{c} $/Å z $ {R}_{c} $/Å z $ {R}_{c} $/Å z I 小腔体 3.875 20 6.593 20 8.056 50 0.973 35.345 大腔体 4.152 21 7.078 24 8.255 50 0.828 14.116 II 小腔体 3.870 20 6.697 20 8.079 20 0.973 35.335 大腔体 4.703 26 7.464 28 8.782 50 2.313 782.847 1.1.2 液相中水的化学电位差
对于纯水相(液态水或冰),Marshall等人提出了
$ \mathrm{\Delta }{\mu }^{W} $ 的计算公式[23]:$$ \frac{\Delta\mu^{\mathrm{\mathit{W}}}}{RT}=\frac{\Delta\mu_{\mathrm{\mathit{W}}}^0}{RT_0}-\int_{T_0}^T\frac{\Delta h_{\mathrm{\mathit{W}}}}{RT^2}\mathrm{d}T+\int_{T_0}^T\frac{\Delta V_{\mathrm{\mathit{W}}}}{RT}\left(\frac{\mathrm{d}p}{\mathrm{d}T}\right)\mathrm{d}T $$ (11) 式中,
$ \mathrm{\Delta }{h}_{W} $ 表示空水合物晶格与纯水相之间的摩尔焓差;$ \mathrm{\Delta }{V}_{W} $ 表示空水合物晶格与纯水相之间的摩尔体积差;$ \mathrm{\Delta }{\mu }_{W}^{0} $ 表示在$ {T}_{0} $ 和$ {p}_{0} $ 条件下空水合物晶格与冰之间的化学位差。对于含烃类溶质富水液相,假定
$ \mathrm{\Delta }{V}_{W} $ 与温度无关,对公式(11)进行了简化[24],如下:$$ \frac{\Delta\mu^{\mathrm{\mathit{W}}}}{RT}=\frac{\Delta\mu_{\mathrm{\mathit{W}}}^0}{RT_0}-\int_{T_0}^T\frac{\Delta h_{\mathrm{\mathit{W}}}}{RT^2}\mathrm{d}T+\int_0^p\frac{\Delta V_{\mathrm{\mathit{W}}}}{RT}\mathrm{d}p-\ln\left(a_{\mathrm{\mathit{W}}}\right)\mathrm{ } $$ (12) 空水合物晶格与纯水相之间的摩尔焓差
$ \mathrm{\Delta }{h}_{W} $ 可以表示为:$$ \begin{aligned} & \Delta h_W=\Delta h_{\mathrm{\mathit{w}}}^0+\int_{T_0}^T\Delta C_{p\mathrm{\mathit{W}}}\mathrm{~d}T \\ & \Delta C_{p\mathrm{\mathit{W}}}=\Delta C_{p\mathrm{\mathit{W}}}^0+b\left(T-T_0\right)\end{aligned} $$ (13) 式中,
$ \mathrm{\Delta }{h}_{W}^{0} $ 为$ {T}_{0} $ 时水在完全空的水合物晶格与纯水相之间的摩尔焓差;$ \mathrm{\Delta }{C}_{pW}^{0} $ 为$ {T}_{0} $ 时水在完全空的水合物晶格与纯水相之间的比热容差;b为比热容的温度系数。$ \mathrm{\Delta }{\mu }_{W}^{0} $ 、$ \mathrm{\Delta }{h}_{W}^{0} $ 、$ \mathrm{\Delta }{C}_{pW}^{0} $ 、$ \mathrm{\Delta }{V}_{W} $ 和b均须通过实验数据回归得到,对不同的水合物结构,取不同的值,如表4所示。表 4 水合物热力学参数参数 I型结构 II型结构 $ \mathrm{\Delta }{\mu }_{W}^{0} $ /(J/mol) 1120 931 $ \mathrm{\Delta }{h}_{W}^{0} $ /(J/mol) 1714 (T<T0) 1400(T<T0) −4297(T>T0) −4611(T>T0) $ \mathrm{\Delta }{V}_{W} $ /(mL/mol) 2.9959 (T<T0) 3.396 44(T<T0) 4.5959(T>T0) 4.996 44(T>T0) $ \mathrm{\Delta }{C}_{pW} $ /(J/mol·K) $ 3.315+0.012(T-{T}_{0}) $(T<T0) $ 1.029+0.004(T-{T}_{0}) $(T<T0) $ -34.583+0.189(T-{T}_{0}) $ (T>T0) $ -36.861+0.181(T-{T}_{0}) $ (T>T0) 1.2 活度系数模型
电解质对气体水合物相平衡的影响主要体现在对富水相中水活度的改变,在不含电解质体系中水的活度近似等于水的摩尔分数。当温度低于冰点时,水的活度
$ {a}_{w}=1 $ ;温度高于冰点时,考虑到溶解在水中的气体和电解质对水活度的影响,水的活度为:$$ a_{w}=a_{w g} a_{ {wel }} $$ (14) 式中,
$ {a}_{wg} $ 表示溶解气对水活度的贡献,awel表示电解质对水活度的贡献。$ {a}_{wg} $ 表达式为:$$ a_{{wg}}=f_{g} \exp \left(A_{g}+B_{g} / T\right) \exp \left[-\frac{\bar{V}_{g}(p-1)}{82.06 T}\right] $$ (15) 式中,
$ {x}_{g} $ 为气体在富水液相中的溶解度;$ \bar{{V}_{g}} $ 为气体分子在富水液相中的偏摩尔体积,C2H4取60,其他取32;$ {A}_{g} $ 和$ {B}_{g} $ 分别为常数,CH4分别取−15.8262、1559.0631。awel表达式为[25]:$$\left\{\begin{array}{l}\ln a_{w e l}=\sum \eta_i m_{e l i} \ln a_{w e l i} / \sum \eta_i m_{e l i} \\ a_{w e l i}=x_w \gamma_w\end{array}\right.$$ (16) 式中,
$ {\eta }_{i} $ 表示第i种电解质中离子的化学计量数。meli表示第i种电解质的摩尔浓度,aweli表示在仅含有第i种电解质的溶液中水的活度,xw表示水的摩尔分数,$ {\gamma }_{w} $ 表示水的活度系数。对于混合电解质,考虑电解质间的相互作用,引入二元相互作用参数,水的活度表示如下:
$$\left\{\begin{array}{l}a_{w e l i}=x_w \gamma_w \xi_{i j} \\ \xi_{i j}=n_1+n_2 \times T+n_3 \times x_{e l}\end{array}\right.$$ (17) 式中,n1、n2和n3是表征电解质之间相互作用的优化参数,见表5。
表 5 表征电解质之间相互作用的优化参数溶液类型 n1 n2 /K−1 n3/(1/摩尔分数) NaCl+LiCl 0.813 6.552e-4 8.725e-3 NaCl+KCl 0.676 2.263e-4 4.962e-3 NaCl+MgCl2 0.729 7.607e-4 6.569e-3 NaCl+CaCl2 0.557 2.673e-4 3.230e-3 NaCl+AlCl3 0.842 5.362e-4 4.236e-3 活度系数是计算水活度的关键。考虑长程相互作用,含盐体系中水的活度系数为[26]:
$$ \ln \gamma=\ln \gamma^{S R I}+\ln \gamma^{L R I} $$ (18) 采用N-NRTL-NRF模型计算短程相互作用[27]:
$$\ln \gamma^{S R I}=x_{e l}^2\left(\zeta_{e l w} \lambda_{e l w}^2+\frac{\zeta_{w e l} \lambda_{w e l}^2}{\beta_{w e l}}-\zeta_{{elw}}-\zeta_{{wel}}\right)$$ (19) 式中,
$ {\zeta }_{elw} $ 和$ {\zeta }_{wel} $ 为模型二元参数,由实验数据计算获取。LiCl、NaCl、KCl、MgCl2、CaCl2和A1Cl3的优化二元参数见表6。表 6 电解质二元优化参数溶液类型 $ {\zeta }_{elw} $ $ {\zeta }_{wel} $ LiCl 11.509 6.375 NaCl 6.359 4.277 KCl 9.620 −10.496 MgCl2 10.629 −3.264 CaCl2 9.259 4.595 AlCl3 18.397 −5.567 $$\begin{cases}\lambda_{w e l}=\dfrac{x_w \beta_{w e l}}{x_w \beta_{w e l}+x_{e l}} & \beta_{w e l}=\exp \left(-\alpha \zeta_{w e l}\right) \\ \lambda_{e l w}=\dfrac{x_{e l} \beta_{e l w}}{x_{e l} \beta_{e l w}+x_w} & \beta_{e l w}=\exp \left(-\alpha \zeta_{e l w}\right)\end{cases}$$ (20) 考虑电解质在液相中的长程相互作用,
$ {\gamma }^{LRI} $ 计算方法为:$$\begin{split} \ln \gamma^{L R I} = & -\left(\frac{1000}{M_w}\right)^{1 / 2} A_\phi\left[\left(\frac{2\left|{\mathrm{z}}_s {\mathrm{z}}_c\right|}{\rho}\right) \ln \left(1+\rho I^{1 / 2}\right)\right. \\& \left. +\left(\frac{\left|{\mathrm{z}}_s {\mathrm{z}}_c\right| I^{1 / 2}-2 I^{3 / 2}}{I+\rho I^{1 / 2}}\right)\right] \end{split}$$ (21) 式中,Mw是水的摩尔质量,
$ {A}_{\phi } $ 是Debye-Hückel常数,$ \rho $ 是最接近参数,I是离子强度。离子强度与离子的电荷数有关,计算公式如下:$$ t=\frac{1}{2}\left(x_{s} {\textit{z}}_{s}+x_{c} {\textit{z}}_{c}\right) $$ (22) 式中,
$ {x}_{c} $ 和$ {x}_{s} $ 为阴、阳离子摩尔分数,$ {{\textit{z}}}_{c} $ 和$ {{\textit{z}}}_{s} $ 为阴、阳离子电荷数。1.3 状态方程
由于多参数状态方程可在更宽T、p范围内准确描述不同物系p-V-T关系,以BWRS状态方程为基础的气-液平衡模型被认为是当前烃类分离计算中最佳的模型之一[28]。BWRS状态方程是一个多参数状态方程,其形式为:
$$\begin{aligned} p =& \rho R T+\left(B R T-A-\frac{C}{T^2}+\frac{D}{T^3}-\frac{E}{T^4}\right) \rho^2 \\ & +\left(b R T-a-\frac{d}{T}\right) \rho^3+\alpha\left(a+\frac{d}{T}\right) \rho^6 \\ & +\frac{c \rho^3}{T^2}\left(1+\gamma \rho^2\right) \exp \left(-\gamma \rho^2\right)\end{aligned}$$ (23) 式中,p表示系统温度,Pa;T表示系统温度,K;ρ表示气相或者液相密度,mol·m−3;R表示气体常数,8.314 J·mol−1·K−1。式中,A、B、C、D、E、a、b、c、d、α、γ为状态方程的11个参数,均可以由它的临界参数Tc、pc、ρc和偏心因子ω从下列公式求得:
$$\left\{\begin{array}{lll}\rho_c B=A_1+B_1 \omega &\qquad \dfrac{\rho_c A}{R T_c}=A_2+B_2 \omega &\qquad \dfrac{\rho_c C}{R T_c^3}=A_3+B_3 \omega \\ \rho_c^2 \gamma=A_4+B_4 \omega &\qquad \rho_c^2 b=A_5+B_5 \omega &\qquad \dfrac{\rho_c^2 a}{R T_c}=A_6+B_6 \omega \\ \rho_c^3 \alpha=A_7+B_7 \omega &\qquad \dfrac{\rho_c^2 c}{R T_c^3}=A_8+B_8 \omega &\qquad \dfrac{\rho_c D}{R T_c^4}=A_9+B_9 \omega \\ \dfrac{\rho_c^2 d}{R T_c^2}=A_{10}+B_{10} \omega &\qquad \dfrac{\rho_c E}{R T_c^5}=A_{11}+B_{11} \omega \exp (-3.8 \omega)\end{array}\right.$$ (24) 利用BWRS状态方程求解逸度系数的表达式如下:
$$\begin{split} \ln \varphi=\ln \frac{f}{p}&= Z-1-\ln Z+\frac{1}{R T} \int_0^p(p-R T \rho) \frac{\mathrm{d} \rho}{\rho^2} \\ & = -\ln \frac{p(V-b)}{R T} + \frac{2}{R}\left(B R-\frac{A}{T}-\frac{C}{T^3}+\frac{D}{T^4}-\frac{E}{T^5}\right) \rho \\ & +\frac{3}{2 R}\left(b R-\frac{a}{T}-\frac{d}{T^2}\right) \rho^2+\frac{6 \alpha}{5 R T}\left(a+\frac{d}{T}\right) \rho^5 \\& + \frac{c}{\gamma R T^3}\left[1+\left(\frac{\gamma \rho^2}{2}+\gamma^2 \rho^4-1\right) \exp \left(-\gamma \rho^2\right)\right]\end{split}$$ (25) 利用上述方程求解气体逸度时,只需要求得气体的摩尔密度即可。气相摩尔密度的求取通常以理想气体密度p/RT为初始值,运用Newton 迭代法进行求解:
$$ \rho_{k+1}=\rho_{k}-\frac{f\left(\rho_{k}\right)}{f^{\prime}\left(\rho_{k}\right)} $$ (26) 为了方便迭代计算,将BWRS状态方程表示如下:
$$ \begin{aligned} f(\rho)= & \rho R T+\left(B R T-A-\frac{C}{T^{2}}+\frac{D}{T^{3}}-\frac{E}{T^{4}}\right) \rho^{2}+ \\ & \left(b R T-a-\frac{d}{T}\right) \rho^{3}+\alpha\left(a+\frac{d}{T}\right) \rho^{6}+ \\ & \frac{c \rho^{3}}{T^{2}}\left(1+\gamma \rho^{2}\right) \exp \left(-\gamma \rho^{2}\right)-p \end{aligned} $$ (27) BWRS状态方程中的相关参数以及各种气体的物理参数如表7和表8所示。
表 7 BWRS状态方程相关参数参数下标j 参数值 Aj Bj 1 0.443 690 0.115 449 2 1.284 380 −0.920 731 3 0.356 306 1.708 710 4 0.544 979 −0.270 896 5 0.528 629 0.349 621 6 0.484 011 0.754 130 7 0.070 523 −0.044 448 8 0.504 087 1.322 450 9 0.030 745 0.179 433 10 0.073 283 0.463 492 11 0.006 450 −0.022 143 表 8 气体的临界参数气体类型 Tc/K pc/MPa ρc/(kmol·m−3) ω N2 126.15 3.394 11.0990 0.0350 CO2 304.09 7.376 10.6380 0.2100 H2O 647.30 22.048 17.8570 0.3440 CH4 190.69 4.604 10.0500 0.0130 C2H6 305.38 4.880 6.7566 0.1018 C3H8 369.89 4.250 4.9994 0.1570 2. 模型求解与验证
2.1 模型求解
气体水合物相平衡条件是依据水在水合物相与在富水相中的化学位相等来确定的。计算水在水合物相的化学位
$ \mathrm{\Delta }{\mu }^{H} $ 有2个关键的参数需要求取:①气相组分的逸度(采用BWRS状态方程);②气体的Langmuir常数(采用Kihara势能模型)。确定水在富水相中的化学位的关键是需要计算在电解质存在条件下水的活度(采用Pitzer-Debye-Hückel模型和N-NRTL-NRF模型)。详细的计算流程如图1所示。采用温度绝对平均相对偏差(AARD-T)评估实验值与预测值的偏差。
$$ A A R D-T=\frac{1}{N} \sum_{i=1}^{N}\left|\frac{T_{i}^{K}-T_{i}^{C}}{T_{i}^{E}}\right| $$ (28) 2.2 模型验证
为了验证模型准确性,收集了250个纯水中甲烷水合物相平衡数据点,并将其与模型预测结果进行对比[12-13,28-29]。同时,分别使用Haghighi[14]、Hsieh[15]和Li等人[17]提出的模型进一步对比分析模型适用范围。如图2所示,温压条件分别在273.2~303.48 K和2.63~72.26 MPa范围内,模型与实验数据吻合程度高,偏差小,其AARD-T仅为0.08。相同压力下,Hsieh模型预测温度大都高于实验值,而Haghighi模型的预测温度则始终低于实验值。且随着压力的升高,Haghighi模型与Li模型的预测温度逐渐偏离实验值,这主要是因为计算Langmuir吸附常数时使用的二元相互作用经验参数不再适用。Haghighi和Li的AARD-T分别为0.20、0.68、0.13。总而言之,所建立的模型能够准确预测纯水条件下甲烷水合物相平衡条件。
收集CH4水合物在不同浓度LiCl、NaCl、KCl、MgCl2、CaCl2和AlCl3溶液中相平衡实验数据,评估模型准确性。如图3(a)-3(f)所示,实验数据与预测的温度-压力-盐浓度三维剖面很好的贴合,预测结果与实验值之间具有很好的一致性,模型的AARD-T比其他模型都要低,如表9所示。
为了进一步分析该模型对天然气水合物在混合电解质溶液中相平衡条件预测的适用性,采用在NaCl+KCl和NaCl+CaCl2混合盐溶液中CH4水合物测试数据对模型进行验证。如图4(a)-4(b)所示,模型预测值与文献中的实验数据基本一致,模型的AARD-T比其他模型都要低,分别为0.04和0.06,表明考虑电解质之间的相互作用是非常必要的。同时也表明该模型可以在相对广泛的温压范围内准确预测盐溶液中天然气水合物的相平衡条件。
表 9 电解质溶液中CH4水合物相平衡预测结果液体类型 浓度范围 T/K P/MPa AARD-T /% N Haghighi Hsieh Li Chin 本文 纯水 273.30~303.48 2.680~72.260 0.68 0.20 0.13 0.35 0.08 250 LiCl 0.050~0.200 (wt) 260.17~288.81 3.506~22.155 0.22 0.42 0.16 0.23 0.15 17 NaCl 0.040~0.090 (me) 261.85~278.05 2.390~11.000 0.26 0.36 0.15 0.30 0.13 23 KCl 0.050~0.100 (wt) 269.16~283.20 1.830~8.820 0.34 0.16 0.12 0.28 0.06 17 MgCl2 0.030~0.150 (wt) 270.75~286.40 2.820~12.955 0.09 0.16 0.15 0.13 0.09 22 CaCl2 0.008~0.053 (me) 260.00~284.40 4.920~10.290 0.25 0.17 0.20 0.62 0.10 26 AlCl3 0.090~0.150 (wt) 272.15~278.15 4.040~8.380 0.21 0.26 0.24 0.34 0.03 8 NaCl+KCl 0.030~0.150 (wt) 263.58~281.46 2.569~9.379 0.11 0.09 0.26 0.52 0.04 37 NaCl+CaCl2 0.030~0.100 (wt) 266.02~281.76 2.504~9.664 0.13 0.06 0.17 0.83 0.06 24 注:wt为质量分数,me为摩尔分数。 3. 结果与讨论
基于温度-压力-盐浓度三维相平衡曲面,分析化学因素对三维相平衡曲面的影响以及相平衡曲面的Clausius-Clapeyron线性行为。
3.1 化学因素对三维相平衡面的影响
图5为不同盐浓度下CH4水合物相平衡的p-T曲线。从图中可以看出,在较高压力下,曲线表现出相似平移关系;但在较低压力条件下,平移关系则不明显。图6为不同电解质下CH4水合物相平衡温度梯度分布示意图。从图6(a-c)中可以看出,在0~20 MPa压力范围内,随着盐浓度的增加,温度梯度波动较为剧烈。表明在相同压力下,CH4水合物相平衡温度随盐浓度的增加变化有显著差异,p-T曲线缺乏平移特征。同时,MgCl2、CaCl2和AlCl3盐溶液,CH4的三维水合物相平衡面非线性平移区扩大,如图6(d-f)所示。从图6还可以看出,在恒压条件下,随着盐浓度的增加,相平衡温度梯度逐渐减小。这表明随着盐浓度的升高,CH4水合物平衡温度降低,更容易发生解离。不同的电解质在相同摩尔分数下表现出不同的抑制效果,总体上6种氯盐对水合物CH4的抑制强度排序为:AlCl3 > MgCl2 > CaCl2 > LiCl > NaCl > KCl,这与Li等人的研究结果一致。值得注意的是,Li、Na和K在元素周期表中属于同一族,Mg和Ca也是如此,它们的排列顺序均是从上到下。AlCl3对水合物CH4的抑制作用更强(如图7所示),主要是由于在相同的摩尔分数下,溶液中较高的有效电荷摩尔分数会导致更大的氢键断裂,从而抑制CH4水合物的形成。
3.2 Clausius-Clapeyron线性特征分析
当不考虑电解质影响时,天然气水合物的相平衡关系表现为Clausius-Clapeyron线性特征[30]。绘制lnp-X-1/T曲面图,结果如图8(a-f)所示。结果表明,在特定的温压范围内,电解质水溶液体系中的水合物相平衡曲线符合Clausius-Clapeyron关系,间接证实了模型的准确性。从图中可以看出,不同电解质条件下的1/T-X-lnp曲面的中部整体都接近一个平面,曲面四周表现出不同程度的非线特征。同时,利用Clausius-Clapeyron关系,可以在适用的温压范围内简化电解质水溶液体系中水合物相平衡的预测,不仅可以减少计算时间,同时也提高了模型可用性。从图8可以看出,与其他电解质相比,NaCl和KCl盐溶液中CH4水合物的三维相平衡面呈现出更大的线性区域。因此,实验数据可用于线性拟合NaCl和KCl盐溶液中较宽温压范围内的CH4水合物相平衡曲线,从而建立简化模型。此外,电解质的抑制性越强,1/T和lnp之间的非线性越强(如图9所示)。
4. 结论
1.基于天然气水合物相平衡热力学理论,结合Pitzer-Debye-Hückel方程和N-NRTL-NRF活度系数模型,考虑电解质之间的相互作用关系,建立了温度-压力-氯盐浓度三维相平衡曲面预测模型,适用的氯盐种类以及浓度和压力区间几乎涵盖了工程中可能遇到的所有情况,尤其是在高压、高浓度区域比现有模型表现更优,适用范围更广。
2.在氯盐摩尔分数大于0.02、压力高于20 MPa时,化学因素将导致p-T曲线近似平移;但压力较小时,相平衡温度沿盐浓度方向的梯度剧烈变化,p-T曲线将不再具有平移特性;AlCl3比LiCl、NaCl、KCl、MgCl2和CaCl2对CH4水合物的抑制作用更强。
3.CH4水合物的1/T-X-lnp三维曲面中部表现出较好的Clausius-Clapeyron线性行为,曲面四周具有一定的非线性特征;不同类型氯盐的1/T-X-lnp三维曲面的非线性程度不一样,氯盐电解质对CH4水合物抑制效果越好,1/T和lnp之间的非线性越强。
-
表 1 气体的Kihara势能参数
气体类型 a /Å $ \sigma $ /Å $ \epsilon /k $ /Å CH4 0.3834 3.1650 154.54 C2H6 0.6760 3.1383 190.80 C3H8 0.8340 3.1440 194.55 N2 0.5290 0.2569 150.03 H2S 0.4920 3.1774 198.53 CO2 0.1773 2.9605 170.97 表 2 水合物分子几何参数[22]
水合物
结构腔型 分子核
间距/Å腔体/
单元壳配位数 H2O分子/
单元壳I 512 3.95 2 20 46 51262 4.33 6 24 II 512 3.91 16 20 136 51264 4.73 8 28 表 3 三层球模型相关参数
结构及腔体类型 第一层 第二层 第三层 $ {n}_{0} $ $ {a}_{0} $ $ {R}_{c} $/Å z $ {R}_{c} $/Å z $ {R}_{c} $/Å z I 小腔体 3.875 20 6.593 20 8.056 50 0.973 35.345 大腔体 4.152 21 7.078 24 8.255 50 0.828 14.116 II 小腔体 3.870 20 6.697 20 8.079 20 0.973 35.335 大腔体 4.703 26 7.464 28 8.782 50 2.313 782.847 表 4 水合物热力学参数
参数 I型结构 II型结构 $ \mathrm{\Delta }{\mu }_{W}^{0} $ /(J/mol) 1120 931 $ \mathrm{\Delta }{h}_{W}^{0} $ /(J/mol) 1714 (T<T0) 1400(T<T0) −4297(T>T0) −4611(T>T0) $ \mathrm{\Delta }{V}_{W} $ /(mL/mol) 2.9959 (T<T0) 3.396 44(T<T0) 4.5959(T>T0) 4.996 44(T>T0) $ \mathrm{\Delta }{C}_{pW} $ /(J/mol·K) $ 3.315+0.012(T-{T}_{0}) $(T<T0) $ 1.029+0.004(T-{T}_{0}) $(T<T0) $ -34.583+0.189(T-{T}_{0}) $ (T>T0) $ -36.861+0.181(T-{T}_{0}) $ (T>T0) 表 5 表征电解质之间相互作用的优化参数
溶液类型 n1 n2 /K−1 n3/(1/摩尔分数) NaCl+LiCl 0.813 6.552e-4 8.725e-3 NaCl+KCl 0.676 2.263e-4 4.962e-3 NaCl+MgCl2 0.729 7.607e-4 6.569e-3 NaCl+CaCl2 0.557 2.673e-4 3.230e-3 NaCl+AlCl3 0.842 5.362e-4 4.236e-3 表 6 电解质二元优化参数
溶液类型 $ {\zeta }_{elw} $ $ {\zeta }_{wel} $ LiCl 11.509 6.375 NaCl 6.359 4.277 KCl 9.620 −10.496 MgCl2 10.629 −3.264 CaCl2 9.259 4.595 AlCl3 18.397 −5.567 表 7 BWRS状态方程相关参数
参数下标j 参数值 Aj Bj 1 0.443 690 0.115 449 2 1.284 380 −0.920 731 3 0.356 306 1.708 710 4 0.544 979 −0.270 896 5 0.528 629 0.349 621 6 0.484 011 0.754 130 7 0.070 523 −0.044 448 8 0.504 087 1.322 450 9 0.030 745 0.179 433 10 0.073 283 0.463 492 11 0.006 450 −0.022 143 表 8 气体的临界参数
气体类型 Tc/K pc/MPa ρc/(kmol·m−3) ω N2 126.15 3.394 11.0990 0.0350 CO2 304.09 7.376 10.6380 0.2100 H2O 647.30 22.048 17.8570 0.3440 CH4 190.69 4.604 10.0500 0.0130 C2H6 305.38 4.880 6.7566 0.1018 C3H8 369.89 4.250 4.9994 0.1570 表 9 电解质溶液中CH4水合物相平衡预测结果
液体类型 浓度范围 T/K P/MPa AARD-T /% N Haghighi Hsieh Li Chin 本文 纯水 273.30~303.48 2.680~72.260 0.68 0.20 0.13 0.35 0.08 250 LiCl 0.050~0.200 (wt) 260.17~288.81 3.506~22.155 0.22 0.42 0.16 0.23 0.15 17 NaCl 0.040~0.090 (me) 261.85~278.05 2.390~11.000 0.26 0.36 0.15 0.30 0.13 23 KCl 0.050~0.100 (wt) 269.16~283.20 1.830~8.820 0.34 0.16 0.12 0.28 0.06 17 MgCl2 0.030~0.150 (wt) 270.75~286.40 2.820~12.955 0.09 0.16 0.15 0.13 0.09 22 CaCl2 0.008~0.053 (me) 260.00~284.40 4.920~10.290 0.25 0.17 0.20 0.62 0.10 26 AlCl3 0.090~0.150 (wt) 272.15~278.15 4.040~8.380 0.21 0.26 0.24 0.34 0.03 8 NaCl+KCl 0.030~0.150 (wt) 263.58~281.46 2.569~9.379 0.11 0.09 0.26 0.52 0.04 37 NaCl+CaCl2 0.030~0.100 (wt) 266.02~281.76 2.504~9.664 0.13 0.06 0.17 0.83 0.06 24 注:wt为质量分数,me为摩尔分数。 -
[1] 刘昌岭, 孙运宝. 海洋天然气水合物储层特性及其资源量评价方法[J]. 海洋地质与第四纪地质,2021,41(5):44-57.LIU Changling, SUN Yunbao. Characteristics of Marine gas hydrate reservoir and its resource evaluation methods[J]. Marine Geology & Quaternary Geology, 2021, 41(5):44-57. [2] 庞雄奇, 胡涛, 蒲庭玉, 等. 中国南海天然气水合物资源产业化发展面临的风险与挑战[J]. 石油学报,2024,45(7):1044-1060. doi: 10.7623/syxb202407002PANG Xiongqi, HU Tao, PU Tingyu, et al. Risks and challenges of the industrial development of methane hydrate resources in the South China Sea[J]. Acta Petrolei Sinica, 2024, 45(7):1044-1060. doi: 10.7623/syxb202407002 [3] 于倩男, 唐慧敏, 李承龙, 等. 天然气水合物注热分解渗流特征及数值模拟[J]. 东北石油大学学报,2023,47(6):38-54. doi: 10.3969/j.issn.2095-4107.2023.06.004YU Qiannan, TANG Huimin, LI Chenglong, et al. Numerical simulation of seepage flow and decomposition characteristics of natural gas hydrate in process of heat injection production[J]. Journal of Northeast Petroleum University, 2023, 47(6):38-54. doi: 10.3969/j.issn.2095-4107.2023.06.004 [4] 邵亚洲. 多孔介质中天然气水合物降压分解实验与数值模拟研究[D]. 哈尔滨: 哈尔滨工程大学, 2022.SHAO Yazhou. Experimental and numerical research ondissociation characteristics of natural gashydrate by depressurization in porous media[D]. Harbin: Harbin Engineering University, 2022. [5] 陈兵兵. 多孔介质内水合物相变与流体运移相互作用机制研究[D]. 大连: 大连理工大学, 2021.CHEN Bingbing. Study on interaction mechanism between fluid migration and hydrate phase transition in porous media[D]. Dalian: Dalian University of Technology, 2021. [6] 王俊杰, 蒋彬, 廖玉宏, 等. 南海琼东南盆地两种不同成因天然气水合物赋存的深层沉积物地球化学特征对比[J]. 地球化学,2023,52(2):180-190.WANG Junjie, JIANG Bin, LIAO Yuhong, et al. Comparison of geochemical characteristics of natural gas hydrate-related sediments in Qiongdongnan Basin, South China Sea[J]. Geochimica, 2023, 52(2):180-190. [7] 李明川, 樊栓狮, 徐赋海, 等. 天然气水合物热分解Stefan相变数学模拟研究[J]. 化工学报,2021,72(6):3252-3260.LI Mingchuan, FAN Shuanshi, XU Fuhai, et al. Mathematical modeling of Stefan phase change for thermal dissociation of natural gas hydrate[J]. CIESC Journal, 2021, 72(6):3252-3260. [8] 吴艳辉, 代锐, 张磊, 等. 深水井筒海水聚合物钻井液水合物生成抑制与堵塞物处理方法[J]. 钻井液与完井液,2023,40(4):415-422. doi: 10.12358/j.issn.1001-5620.2023.04.001WU Yanhui, DAI Rui, ZHANG Lei, et al. Study on methods of gas hydrate inhibition and blockage treatment using seawater polymer muds in deepwater well drilling[J]. Drilling Fluid & Completion Fluid, 2023, 40(4):415-422. doi: 10.12358/j.issn.1001-5620.2023.04.001 [9] 任冠龙, 孟文波, 何玉发, 等. 深水浅层钻井液水合物抑制性能优化[J]. 钻井液与完井液,2022,39(5):529-537. doi: 10.12358/j.issn.1001-5620.2022.05.001REN Guanlong, MENG Wenbo, HE Yufa, et al. Optimization of hydrate inhibition performance of deep water shallow drilling fluid[J]. Drilling Fluid & Completion Fluid, 2022, 39(5):529-537. doi: 10.12358/j.issn.1001-5620.2022.05.001 [10] 张更. 海洋浅表层天然气水合物固态流化开采井筒多相流动规律研究[D]. 北京: 中国石油大学(北京), 2023.ZHANG Geng. Research on the wellbore multiphase flow rules during natural gas hydrate solid fluidization exploitation in offshore shallow surface[D]. Beijing: China University of Petroleum(Beijing), 2023. [11] JAVANMARDI J, MOSHFEGHIAN M, MADDOX N R. An accurate model for prediction of gas hydrate formation conditions in mixtures of aqueous electrolyte solutions and alcohol[J]. The Canadian Journal of Chemical Engineering, 2001, 79(3):367-373. doi: 10.1002/cjce.5450790309 [12] MOHAMMADI H A, ANDERSON R, TOHIDI B. Carbon monoxide clathrate hydrates: Equilibrium data and thermodynamic modeling[J]. AIChE Journal, 2005, 51(10):2825-2833. doi: 10.1002/aic.10526 [13] MOHAMMADI A H, TOHIDI B. Prediction of hydrate phase equilibria in aqueous solutions of salt and organic inhibitor using a combined equation of state and activity coefficient‐based model[J]. The Canadian Journal of Chemical Engineering, 2005, 83(5):865-871. [14] BANDYOPADHYAY A A, KLAUDA J B. Gas hydrate structure and pressure predictions based on an updated Fugacity-Based model with the PSRK equation of state[J]. Industrial & Engineering Chemistry Research, 2011, 50(1):148-157. [15] HAGHIGHI H, CHAPOY A, TOHIDI B. Methane and water phase equilibria in the presence of single and mixed electrolyte solutions using the Cubic-Plus-Association equation of state[J]. Oil & Gas Science and Technology, 2009, 64(2):141-154. [16] HSIEH M K, YEH Y T, CHEN Y P, et al. Predictive method for the change in equilibrium conditions of gas hydrates with addition of inhibitors and electrolytes[J]. Industrial & Engineering Chemistry Research Journal, 2012, 51(5):2456-2469. [17] MA Q L, CHEN G J, GUO T M. Modelling the gas hydrate formation of inhibitor containing systems[J]. Fluid Phase Equilibria, 2003, 205(2):291-302. doi: 10.1016/S0378-3812(02)00295-9 [18] LI S G, LI Y J, WANG J Q, et al. Prediction of gas hydrate formation conditions in the presence of electrolytes using an N-NRTL-NRF activity coefficient model[J]. Industrial & Engineering Chemistry Research, 2020, 59(13):6269-6278. [19] Parrish W, PARRISH J. Dissociation pressures of gas hydrates formed by gas mixtures[J]. Industrial & Engineering Chemistry Research, 1972, 11(1):26-35. [20] JOHN V T, PAPADOPOULOS K D, HOLDER G D. A generalized model for predicting equilibrium conditions for gas hydrates[J]. AIChE Journal, 2010, 31(2):252-259. [21] MCKOY V, SINANOĞLU O. Theory of dissociation pressures of some gas hydrates[J]. The Journal of Chemical Physics, 1963, 38(12):2946-2956. doi: 10.1063/1.1733625 [22] MEHTA A P, SLOAN E D. Improved thermodynamic parameters for prediction of structure H hydrate equilibria[J]. AIChE Journal, 1996, 42(7):2036-2046. doi: 10.1002/aic.690420724 [23] SAITO S, MARSHALL D R, KOBAYASHI R. Hydrates at high pressures: Part II. Application of statistical mechanics to the study of the hydrates of methane, argon, and nitrogen[J]. AIChE Journal, 1964, 10(5):734-740. doi: 10.1002/aic.690100530 [24] HOLDER G D, CORBIN G, PAPADOPOULOS K D. Thermodynamic and molecular properties of gas hydrates from mixtures containing methane, argon, and krypton[J]. Industrial & Engineering Chemistry Fundamentals, 1980, 19(3):282-286. [25] NASRIFAR K, MOSHFEGHIAN M. A model for prediction of gas hydrate formation conditions in aqueous solutions containing electrolytes and/or alcohol[J]. The Journal of Chemical Thermodynamics, 2001, 33(9):999-1014. doi: 10.1006/jcht.2000.0811 [26] HSIEH C M, SANDLER S I, LIN S T. Improvements of COSMO-SAC for vapor–liquid and liquid–liquid equilibrium predictions[J]. Fluid Phase Equilibria, 2010, 297(1):90-97. doi: 10.1016/j.fluid.2010.06.011 [27] HAGHTALAB A, SHOJAEIAN A, MAZLOUMI S H. Nonelectrolyte NRTL-NRF model to study thermodynamics of strong and weak electrolyte solutions[J]. The Journal of Chemical Thermodynamics, 2011, 43(3):354-363. doi: 10.1016/j.jct.2010.10.004 [28] BÖTTGER A, KAMPS Á P S, MAURER G. An experimental investigation of the phase equilibrium of the binary system (methane plus water) at low temperatures: solubility of methane in water and three-phase (vapour plus liquid plus hydrate) equilibrium[J]. Fluid Phase Equilibria, 2016, 407:209-216. doi: 10.1016/j.fluid.2015.03.041 [29] MAEKAWA T. Equilibrium conditions for clathrate hydrates formed from methane and aqueous propanol solutions[J]. Fluid Phase Equilibria, 2008, 267(1):1-5. doi: 10.1016/j.fluid.2008.02.006 [30] 孙始财, 刘昌岭, 业渝光, 等. 氯盐溶液中甲烷水合物高压分解条件及影响因素[J]. 物理化学学报,2011,27(12):2773-2778. doi: 10.3866/PKU.WHXB20112773SUN Shicai, LIU Changling, YE Yuguang, et al. Dissociation conditions and influencing factors of methane hydrate in chloride salt solution under high pressure[J]. Acta Physico-Chimica Sinica, 2011, 27(12):2773-2778. doi: 10.3866/PKU.WHXB20112773 -