不同相体系下H2S水合物生长过程分子动力学模拟

罗泽东，李爱蓉，刘发平，陆 原

（1.西南石油大学 化学化工学院，四川 成都 610500；2.中海油（天津）油田化工有限公司，天津 300450）

目前我国已探明气田中有一百多个高含硫气田（H2S体积分数在2%～10%之间），其累计探明储量已超过1012m3，其中普光气田天然气中H2S含量高达13%～18%[1]。高含硫天然气比普通天然气在较高温度和较低压强下更容易形成气体水合物[2]。其原因是在天然气中H2S气体的存在会促进水合物的形成，且一旦形成水合物就很难去除，在生产、运输和加工过程中造成管道堵塞、采气管道腐蚀和安全隐患等问题。因此，研究H2S水合物的形成过程尤为重要。但是，水合物的成核与生长过程是一个微观现象，很难通过实验进行观察。同时，H2S具有腐蚀性和剧毒性使得从实验层面研究其水合物的形成过程更具有挑战性。

目前，分子动力学模拟（MD模拟）可以从分子层面考察气体水合物成核与生长的微观规律[3-8]。近年来已有文献报道借助MD模拟探究了纯H2S水合物的成核与生长行为。Liang等[9]模拟了在230～250 K的温度范围内H2S水合物的微观成核过程，结果表明，相对于CH4水合物，H2S水合物的成核速率更大，其原因在于H2S气体在水中具有更大的溶解度。Liang等[10]采用局部过热方法，熔化一半H2S水合物晶体，进而随机产生H2S气相、H2S溶液相和H2S水合物固相的三相体系与H2S溶液相、H2S水合物固相的两相体系；再采用局部过热与局部过冷的温度调控方式，分别研究两个体系下H2S水合物的生长过程。模拟结果表明，三相体系和两相体系下H2S水合物的生长速率远远大于CH4水合物，这主要归因于H2S分子在水中的高溶解度和迁移率。但是，从初始模型与研究手段来看，这种局部过热耦合过冷的方式明显加速了H2S水合物的生长速率。同时，水合物晶体种子也会促进水合物笼对溶液中H2S的吸附。这可能掩盖了不同体系下H2S水合物更为真实的生长过程，所揭示的H2S水合物生长的微观机理并不完善。

为了进一步探究温度对不同体系下H2S水合物生长过程的影响，本文通过MD模拟，分别在H2S气相、溶液相和水合物相的三相体系，以及H2S溶液相和水合物相的两相体系中，研究了H2S水合物在较宽温度范围内的晶体生长过程。就生长厚度、扩散系数、局部密度等因素对H2S水合物生长速率的影响进行了深入分析。

1 模型建立与模拟方法

1.1 模型建立

1.1.1 气液固三相模型

H2S水合物属于Ⅰ型水合物，晶格常数为1.1877 nm。本文H2S水合物单晶胞的氧原子坐标取自X射线衍射实验结果[11]，然后对单晶胞的3个方向扩胞形成2×2×10单元的超晶体，并将其放置于2.4 nm×2.4 nm×12.0 nm的正交模拟盒子（z方向长度为12.0 nm）。本文在Liang等[7]关于H2S水合物生长模型的基础上，只保留少部分晶体种子，以此建立三相与两相体系下H2S水合物的生长模型。其中，三相建模细节为：（1）固定z方向左边（2.4 nm × 2.4 nm × 2.4 nm）三维区域中的水合物晶体，剩余部分在350 K和10 MPa，等温等压系综（NPT系综）中，熔化2.5 ns后，产生H2S气相、H2S溶液相和H2S水合物固相三相体系；（2）继续固定左边的晶体，熔化区从350 K分别降温至285 K、275 K、270 K、265 K、260 K和250 K，降温速率为0.75 K/ps，使熔化区逐渐平衡至各生长温度；（3）随后固定熔化区，左边晶体分别在各个温度下恒温100 ps，以释放晶体区与熔化区之间的应力。

形成的H2S气相、溶液相、水合物固相的三相初始模型如图1所示。本文所有的MD模拟程序采用LAMMPS[12]软件进行。

图1 H2S气液固三相初始MD模型

1.1.2 液固两相模型

为了使初始构型中不包含H2S气相，故随机移除H2S水合物单晶胞中一半的H2S分子，然后对单晶胞的3个方向扩胞形成2 × 2 × 10单元的H2S水合物超晶体，此时超晶体的晶穴占有率为50%。并将其放置于2.4 nm × 2.4 nm × 12.0 nm的正交模拟盒子（z方向长度为12.0 nm）。熔化和平衡过程的模拟细节与三相模型相同，H2S溶液相与H2S水合物固相的两相初始模型如图2所示。

图2 H2S液固两相初始MD模型

1.2 模拟方法

水分子间的相互作用势选用TIP4P/2005势函数[13]，而H2S分子之间的相互作用采用Kamath[14]三点模型，如表1所示，这些参数能够准确反映H2S水合物的相平衡，并且在以往的研究中得到了广泛的应用[4,5]。H2O分子与H2S分子的Lennard-Jones（LJ）交叉作用参数由Lorentz-Berthelot混合规则计算，长程库伦作用采用pppm算法[15]求解，LJ交叉作用与长程库伦作用的截断半径均选取1.0 nm。在整个模拟过程中应用三维周期边界条件，时间步长为2 fs。对于产出阶段设定NPT系综，在6个初始温度和10 MPa下正式模拟100 ns，分别研究不同体系下H2S水合物的生长情况。在每个温度下均进行3次不同初速度随机数种子的模拟，以求取数据平均值。

表1 LJ势函数参数及原子电荷分配

采用3种不同的措施来跟踪H2S水合物的生长过程，具体为：（1）氧原子的径向分布函数（RDF）与H2S分子的配位数；（2）H2S和H2O分子的均方位移曲线（MSD）；（3）H2S分子密度分布。

RDF是指距离任意参考粒子α中心r处出现β粒子的几率，β粒子可以与α粒子相同也可以不同，可以定义为式(1)。对g(r)在[0,r0]进行积分，可计算得到参考粒子在r0范围内的临近配位数。

式中，dN表示r→r+dr间的粒子数目；ρ为系统的密度，g/cm3。

均方位移如式(2)：

式中，MSD(t)为分子的均方位移，10-20m2；尖括号表示对总原子数的平均；ri(t)和ri(0)分别表示i粒子在t时刻的位置和在0时刻的位置，其中时间t的单位为ns。

对于固体分子来说，运动受到限制，其均方位移在0附近波动；液体分子的运动相对自由，其均方位移随模拟时间增加而增大。根据爱因斯坦的扩散定律，有式(3)：

式中，D为分子扩散系数，10-11m2/s。

1.3 模型验证

图3(a)是温度为260 K、压力为10 MPa下氧原子的RDF。第一波峰在0.1～100 ns内均位于0.273 nm，表明水合物晶胞中五边形的边长为0.273 nm，与文献数据吻合[3,4]。图3(b)为260 K、10 MPa下与H2S分子相关的临近配位数。在距离中心H2S分子的0.572 nm范围内所有的H2O分子被定义为水临近配位数，同时在距离中心H2S分子的0.912 nm范围内所有的H2S分子被定义为硫临近配位数[6]。从图3(b)可以看出，H2S分子的水临近配位数达到20个，表明构成512小笼中水分子个数为20个，与水合物晶体的几何结构保持一致。硫临近配位数逐渐减小至13，H2S分子的第一临近硫配位数与文献[6]保持一致，表明离中心H2S分子的临近H2S分子达到13个。由此可见，在高浓度H2S气相区，H2S气体被消耗用于水合物晶体的生长。

图3 260 K、10 MPa下氧原子径向分布函数(a)与H2S分子的配位数(b)演变

2 结果与讨论

2.1 气液固三相体系下H2S水合物的生长过程

2.1.1 构型分析与生长速率

通过NPT模拟，在285 K、275 K、270 K、265 K、260 K和250 K的初始温度下探究三相体系下H2S水合物晶体生长的微观过程。图4为该体系下H2S水合物100 ns时的生长快照。在100 ns的时间跨度内，6个初始温度在相平衡稳定区以内均发生晶体生长，并且从图4发现在3个较低温度下，最终构型在界面处出现无定型的水笼结构和带有局部特征的H2S气泡；而在图4的3个较高温度下，高温加剧了H2S气泡的聚集程度，最终体系具有明显相分离现象，新晶相呈现长程有序分布。

图4 H2S水合物在三相体系不同初温下100 ns时的生长快照

为了更深入理解三相体系下水合物晶体生长过程中的水分子运动情况，本文在250 K下，随机追踪一个H2S分子在20～80 ns内周围水分子的微观运动现象，如图5所示。

图5 某个H2S分子的周围水结构在20～80 ns的演变过程

由图5可知，在第27 ns时，18个水分子通过氢键作用形成半笼型包围这个H2S分子，36 ns后，周围水分子形成比较完整的512笼型。在48 ns后，512笼继续吸引周围水分子，以此笼为基础连接一个六边形。对于水合物笼的拓扑结构来说，水分子构成的六边形是形成大笼的初始条件。然而，在63 ns时六边形消失，多个水分子形成更复杂的无定型团簇。即使这里的512小笼在10 ns间隔内并未消失，但在远离界面处的液相中，长达数十纳秒并未出现完整的单晶胞水合物。这是因为虽然低温可以提供较大的驱动力促使晶体成核，但三相体系下水分子运动受到限制，扩散较慢，水分子需要频繁更新自身坐标才能到达成核位点，因此需要更长时间生长出单晶胞。该温度下生成的无定型团簇是一个亚稳态中间体，无定型团簇的形成将延迟晶体的进一步生长。这与Jacobson等[16,17]认为高驱动力下，非晶核（即无定型团簇）在动力学上可能比晶核更有利的结果保持一致，因为非晶核形成的自由能势垒较低。

在三相体系下，温度对H2S水合物生长速率的影响表现为其生长速率先增大后减小的趋势，其结果如表2所示。250 K时水合物的生长厚度为3.6 nm，为最小值；260 K时水合物的生长厚度为6.4 nm，为最大值，此时水合物的最大生长速率为0.120 nm/ns，且在260～270 K下，H2S水合物的生长速率相差不大。文献报道的MD模拟与实验结果均表明水合物新晶体的最大生长速率会出现在某个中等温度区间内[18,19]，这说明在三相体系下温度对气体水合物生长过程的影响具有相似性。本文计算得到的H2S水合物生长速率（260 K，0.120 nm/ns）大于同体系下CH4水合物的生长速率（260 K，0.054 nm/ns）[20]。

表2 三相体系下H2S水合物于不同温度的生长数据

2.1.2 扩散系数

为了考察H2S分子和H2O分子在三相体系下运动的快慢，分别计算了H2O和H2S分子在不同温度下的均方位移，结果如图6所示。由图6可知，随生长温度升高，H2O与H2S分子的MSD值逐渐增大。与H2O分子相比，20 ns后H2S分子的MSD曲线出现较大波动，表明高温加速H2S分子运动，使溶液相中的H2S分子朝向H2S气相区运动，相分离状态更加明显，从而延缓晶体的进一步生长。

根据式(3)，通过拟合图6中H2O与H2S分子的MSD曲线，得出各温度下H2O与H2S分子的扩散系数，如表3所示。H2S分子的扩散系数约为H2O分子扩散系数的8～10倍，且H2O分子的扩散系数（275 K，2.472 ×10-10m2/s）小于文献中纯水的自扩散系数（275 K，11.49 × 10-10m2/s）[21]，表明晶体生长过程中水的运动受到限制，H2S水合物的生长速率主要受H2S分子扩散的影响。而文献报道关于CH4水合物的三相生长过程中，发现CH4的扩散系数小于H2O的扩散系数[20,22]。可见，H2S分子在水相中的扩散能力大于CH4分子，使得三相体系下H2S水合物的生长速率大于CH4水合物。因此气体水合物的生长速率与气体分子在水相中的扩散能力较强相关。

图6 三相体系下H2O分子(a)与H2S分子(b)均方位移随时间的变化

表3 三相体系下不同温度H2O分子与H2S分子的扩散系数

2.1.3 H2S的密度分布

通过计算H2S分子在z方向的质量密度分布可以观察出新晶体在空间中排布的对称性，如图7所示。

图7 三相体系下z方向H2S分子的密度分布

由图7可知，在285 K下，新生长出的晶体中H2S密度分布呈现上下波动、周期性复现趋势，与晶体种子区的密度变化趋势保持一致。而250 K下，H2S分子的密度出现局部最大值，表明低温下较低的扩散系数使得密度分布不均。初始温度增加可以强化H2O分子与H2S分子的扩散，使得新晶相呈现长程有序分布，与图4不同温度下的水合物生长快照保持一致。

2.2 液固两相体系下H2S水合物的生长过程

2.2.1 构型分析与生长速率

从图8可以看出，250 K下的最终构型呈现无定型的水笼结构（黄色方框区域）；其余5个温度的最终构型表现为新晶相长程有序，仅有少数H2S分子溶解在水中。与三相体系不同的是，在较高温下生长构型并未出现相分离状态。

图8 H2S水合物在两相体系不同初温下50 ns时的生长快照

表4为不同初始温度下H2S水合物的两相生长数据。250 K下水合物生长厚度为3.0 nm，为最小值；260 K下水合物生长厚度为6.9 nm，达到最大值。除250 K外，其余5个生长温度对应的生长速率相差不大，约为0.220 nm/ns。与三相体系相比，两相体系下温度对H2S水合物生长过程的影响不是很敏感，在初始温度为260 K条件下，两相体系下生长速率0.228 nm/ns大于三相水体系的生长速率0.120 nm/ns。

表4 两相体系下H2S水合物于不同温度的生长数据

为进一步观察两相体系的H2S水合物晶体生长，跟踪了260 K下固液界面处晶体厚度的演变过程。图9为两相体系下0.5～12 ns内H2S水合物晶体生长过程，第7 ns时，新生长出的厚度为1.6 nm，瞬时生长速率可达0.228 nm/ns，该体系下表现为界面处的逐层生长。表明两相体系下的界面迁移行为可促进H2S水合物的生长，而三相体系下带有局部特征的H2S气泡却延缓了晶体生长。

图9 两相体系下H2S水合物在固液界面处晶体厚度的演变

2.2.2 扩散系数

图10为两相体系下H2O与H2S分子均方位移曲线MSD随时间的变化。

图10 两相体系下H2O分子(a)与H2S分子(b)均方位移随时间的变化

图10(a)中H2O分子的MSD随着时间的变化迅速增加；而从图10(b)可以看出，10 ns后，H2S分子的均方位移曲线变得较为平缓，这是因为两相体系下H2S水合物的生长速率较大，溶液中的H2S分子在较短时间内被水合物笼捕捉完成，最终体系变为H2S水合物和纯水两相。

通过拟合图10中H2O与H2S分子的MSD曲线，得出两相体系下H2O和H2S分子的扩散系数，如表5所示。

表5 两相体系下不同温度H2O分子与H2S分子的扩散系数

两相体系下，H2O分子和H2S分子的扩散系数十分接近。表明两相体系中主体分子与客体分子的扩散能达到相互协同的效果，进而促进H2S水合物晶体生长，增加了H2S水合物的生长速率。

2.2.3 H2S的密度分布

图11为两相体系下100 ns时，z方向H2S分子的密度分布。与250 K相比，285 K下z方向H2S分子密度峰值较大且分布有序。这表明，通过增加初始温度，可使得新生长出的水合物晶体变得长程有序，与三相体系下H2S分子的密度分布相似。

图11 两相体系下z方向H2S分子的密度分布

3 结论

本文采用MD模拟，研究了三相和两相体系下温度对H2S水合物生长过程的影响，并借助生长构型、均方位移曲线和密度分布等手段对比讨论了不同温度条件下H2S水合物生长过程的差异，得到主要结论如下：

（1）三相体系下，温度对H2S水合物生长过程的影响较大。无定型结构出现在低温区域，而具有长程有序的新晶体出现在高温区域。250～285 K温度范围内，H2S水合物的最大生长速率会出现在某个中等温度区间内。

（2）两相体系下，温度对H2S水合物生长过程的影响不敏感。除250 K，其余温度下生长出的新晶体均表现为长程有序，且生长速率相差不大，约为0.220 nm/ns。

（3）H2S水合物在两相体系下的生长速率大于三相体系。三相体系下H2S分子的扩散系数约为H2O分子扩散系数的8～10倍，较大的扩散系数使得体系产生相分离状态，不利于晶体生长；而在两相体系下H2O分子与H2S分子的扩散能达到相互协同的效果，促进晶体生长。

本研究结论可为含硫气藏天然气开采中水合物的防治提供理论参考。尤其是H2S分子溶解在水相的体系，可以考虑添加能够减弱H2O与H2S分子间相互扩散能力的抑制剂，在水合物生长初期抑制其形成。