深海已经成为世界油气资源开发的重要接替区域，但深海高压低温条件下形成的天然气水合物将严重影响钻井完井作业及海洋装备的安全^[1–6]。实践表明，使用动力学抑制剂可有效抑制深水钻井过程中天然气水合物的生成，从而降低钻井风险。目前，对抑制剂作用机理的主要认识是，抑制剂分子通过氢键吸附在天然气水合物笼形结构表面，阻止天然气水合物笼有效成长或导致其畸形发展。但是抑制剂作用一段时间以后，仍然会形成大量的天然气水合物，上述认识无法解释这种现象。因此，有必要研究天然气水合物动力学抑制剂与水分子的相互作用情况，进一步明确其相互作用机理。

目前国内外主要通过模拟方式进行相关研究，对天然气水合物动力学抑制剂作用机理的研究大多集中于抑制剂分子官能团、天然气水合物结构微观变化、水分子参与成键等方面：D. M. Duffy等人^[7]发现PVP分子可以通过扰乱水分子的团簇结构，阻止天然气水合物晶核的形成，但是其对已形成水合物晶粒的作用微乎其微；B. Kvamme等人^[8–9]提出用分子模拟条件下体系自由能的变化幅度作为评价动力学抑制剂效果的参考手段，体系自由能越高，则动力学抑制剂的抑制效果越好；M.T.Storr等人^[10]模拟了抑制剂分子在Ⅱ型天然气水合物表面的吸附，发现季铵盐类动力学抑制剂可以通过氢键作用吸附在空腔之上，从而阻止客体分子的进入。笔者在参考前人研究成果的基础上，通过分子模拟，建立了抑制剂溶液–甲烷水合物体系模型（以下称为体系模型），对含有环状结构的PVP和PVCap，含有链状结构的PMC，以及同时含有环状结构和链状结构的新型动力学抑制剂YZ^[11]与水分子之间的相互作用进行了模拟研究，主要包括体系中水分子的扩散系数、抑制剂分子与水分子形成的氢键数目及氢键能的变化规律，以期明确不同分子结构动力学抑制剂与水分子之间的相互作用规律。

1.   体系模型的建立及模拟过程

采用I型甲烷水合物晶胞进行模拟，晶胞中氧原子的位置参照X射线衍射实验结果^[12]。建模时，氧原子依照坐标添加，氢原子按照Bernal-Fowler规则^[13]随机分布在氧原子周围。建立I型天然气水合物3×3×3晶胞后，切割天然气水合物晶胞的（0 0 1）面^[14]作为模型的底部，在天然气水合物晶胞上方添加800个水分子作为模型的中部，模型顶部为真空层。体系模型的构象如图1所示（图1中，红色、白色和灰色的棍棒分别为氢原子、氧原子和碳原子，蓝色的虚线为天然气水合物晶胞的氢键）。

图  1  抑制剂溶液-甲烷水合物体系模型棍棒状构象

Figure  1.  Stick-shaped conformation of inhibitor solution-methane hydrate system model

下载: 全尺寸图片 幻灯片

模拟在Material Studio软件中进行。Material Studio软件是专门为材料科学领域研究开发的一款可在PC机上运行的模拟软件。利用该软件能很好地研究各种晶体、无定型以及高分子材料的性质与相关过程。选择一致性价力场 (consistent valence force field，CVFF)描述天然气水合物分子之间的相互作用^[15–16]。模拟在NVT系综（等粒子数、等体积、等温度环境）下进行，并考虑周期性边界条件，用Lennard-Jones势模型^[17]处理甲烷分子与水分子之间的相互作用，采用Ewald方法^[17]处理长程静电作用，截断半径为1.25 nm；采用Nose-Hoover方法^[17]控制温度，弛豫时间为1 ps，时间步长为0.5 fs，模拟时间为250 ps。

PVP，PVCap，PMC和YZ等4种动力学抑制剂是由不同种类单体聚合而成的聚合物，聚合物的链节是指组成聚合物的每一基本重复结构单元。在进行动力学模拟时，分别在建立的基础模型中添加不同数量含有不同数量链节的动力学抑制剂分子，抑制剂分子1～5个，链节1～6个。模拟过程分为2步：1）固定天然气水合物结构，模拟15 ps，释放天然气水合物与无机盐溶液界面处的额外应力；2）松弛固定层之外的天然气水合物，模拟250 ps。观察天然气水合物与动力学抑制剂分子的相互作用过程，提取前200 ps的相关参数，分析动力学抑制剂分子的作用机理。

为了模拟天然气水合物在不同状态时动力学抑制剂与水分子的相互作用，需要先评价模型中天然气水合物的稳定性，确定该水合物处于不同状态时的温度。自扩散系数是反映体系稳定性的重要参数之一，自扩散系数越低，体系越趋于稳定；反之，自扩散系数越高，体系内各分子的扩散性越强，体系越不稳定。自扩散系数可由均方位移（MSD）曲线的Einstein方程^[18]求得，MSD曲线反映的是动力学模拟过程中体系分子的运动性质，可用其区分体系的相态。

通过动力学模拟不同温度下天然气水合物基础模型的形态，观察模型底部天然气水合物晶胞的分解状况，并综合考虑实际模拟时长等因素，最终确定：天然气水合物基础模型稳定存在的温度为245.2 K，此时晶胞稳定存在，且没有分解现象；天然气水合物的微分解状态温度为253.2 K，此时晶胞上部半笼形结构发生一定程度的坍塌；天然气水合物的不稳定状态温度为258.2 K，此时晶胞上部半笼形结构完全坍塌，甲烷分子逃逸。通过模拟分析可知，该天然气水合物基础模型稳定存在的温度为245.2 K，水合物的微分解状态温度为253.2 K，水合物的不稳定状态温度为258.2 K。这3种温度下溶液的扩散系数分别为0.15，0.19和0.26 nm²/s，模型中对应的天然气水合物的均方位移如图2所示，其扩散系数分别为0.000 94，0.011 16和0.028 51 nm²/s。

图  2  不同温度下基础模型中天然气水合物晶胞分子的均方位移曲线

Figure  2.  Mean square displacement curve of gas hydrate unit cell molecules in the basic model under different temperatures

下载: 全尺寸图片 幻灯片

2.   模拟结果分析

2.1   扩散系数

动力学抑制剂与溶液中的水分子相互作用，必然会导致溶液的扩散能力发生变化。因此，需要研究溶液自扩散系数的变化。

为了分析相同质量分数动力学抑制剂溶液中，各动力学抑制剂对溶液中水分子扩散系数的影响，定义某一抑制剂（聚合物）分子的名称为X_a-b（X代表动力学抑制剂的类别，a代表该动力学抑制剂的链节数量，b代表溶液中该动力学抑制剂的分子数量）。

选择以下3种情况进行分析（括号中的数据为动力学抑制剂在溶液中的质量分数）：

1）当质量分数约为3.00%时，符合条件的情况有PVP_1-4（3.04%）、PVP_2-2（3.02%）、PVP_4-1（3.00%）、PVCap_1-3（2.85%）、PVCap_3-1（2.83%）、PMC_1-5（2.93%）、PMC_5-1（2.88%）、YZ_1-1（3.23%）；

2）当质量分数约为8.50%时，符合条件的情况有PVP_3-4（8.51%）、PVP_4-3（8.50%）、PVP_6-2（8.49%）、PVCap_2-5（8.86%）、PVCap_5-2（8.83%）、PMC_3-5（8.19%）、PMC_4-4（8.68%）、PMC_5-3（8.17%）、YZ_1-3（9.09%）、YZ_3-1（9.07%）；

3）当质量分数约为10.50%时，符合条件的情况有PVP_3-5（10.42%）、PVP_5-3（10.40%）、PVCap_3-4（10.43%）、PVCap_4-3（10.41%）、PVCap_6-2（10.40%）、PMC_4-5（10.61%）、PMC_5-4（10.60%）。

由于天然气水合物动力学抑制剂在使用过程中很少会使用聚合物的单体，因此对于PVP、PVCap和PMC抑制剂，应该排除单体数量为1个的情况。当动力学抑制剂的质量分数分别约为3.00%、8.50%和10.50%时，各动力学抑制剂溶液的自扩散系数如图3所示。

图  3  各动力学抑制剂不同质量分数溶液的自扩散系数

Figure  3.  Self-diffusion coefficients of solutions with different mass fractions for each kinetic inhibitor

下载: 全尺寸图片 幻灯片

由图3可知：动力学抑制剂的质量分数大体相同时，具有环状结构的PVP和PVCap抑制剂溶液的自扩散系数要低于具有链状结构的PMC抑制剂溶液；同时，具有环状结构和链状结构的YZ抑制剂溶液的自扩散系数反而更高，这与动力学抑制剂分子不能在水溶液中有效分散有关。整体而言，研究范围内的同一种动力学抑制剂在近似相同的质量分数下，动力学抑制剂分子的链节数量越少，动力学抑制剂分子越容易均匀分散在水中，其水溶液的自扩散系数也相应越低。

2.2   氢键

通常把与氢原子（H）成键的原子称为给体，用D表示；将与氢形成氢健的原子称为受体，用A表示。氢键给体、受体及变量如图4所示（图4中：r_AH，r_AD分别表示受体与氢、受体与给体间的距离，nm；θ和θ_hb为2种用不同方式表示的氢键键角，二者互为补角，（°））。

图  4  氢键给体、受体及变量示意

Figure  4.  Schematic diagram of hydrogen bond donors, receptors and variables

下载: 全尺寸图片 幻灯片

氢键是一种非化学键，形成氢键的原子间既有范德华力作用，又有静电力作用，没有统一标准界定氢键形成的条件，也无法采用一致的函数形式计算氢键的能量。一般认为氢键的形成条件是^[13]：r_AH<0.25 nm，θ_hb<90°（或θ>90°）。

2.2.1   氢键数目

假设一个动力学抑制剂分子理论上可以与水分子或天然气水合物的微晶颗粒形成m个氢键，实际形成了n个氢键，则定义它们的偏差比为：

式中：p为偏差比。

当偏差比大于0时，偏差比越接近于0，表明动力学抑制剂分子与水分子或天然气水合物的微晶颗粒形成的氢键数目与理论值相差越小；当偏差比小于0时，则意味着不仅动力学抑制剂分子中具有孤对电子的原子参与形成氢键，而且其他原子（例如氢原子）也参与形成氢键。分析偏差比有助于对比分析动力学抑制剂结构对抑制剂参与形成氢键的影响。

利用式（1）计算了动力学抑制剂分子不同温度下的偏差比，结果如图5所示。

图  5  不同温度下不同链长的抑制剂分子与水分子形成氢键的偏差比分布

Figure  5.  Deviation ratio distribution of hydrogen bonds formed by inhibitor molecules of different chain lengths and water molecules under different temperatures

下载: 全尺寸图片 幻灯片

由图5可知：温度一定时，具有环状结构的PVP分子和PVCap分子与溶液中水分子实际形成的氢键数远小于其理论值；具有链状结构的PMC分子实际可以与水分子形成数目更多的氢键；同时具有环状结构和链状结构的YZ分子与溶液中水分子实际形成的氢键数同理论值较为接近。抑制剂分子中的氧原子并非都参与形成氢键，如位于分子结构中部的氧原子通常不能参与形成氢键，而位于分子结构（或支链）末端的氧原子才有可能参与形成氢键。温度升高，抑制剂分子与水分子形成的氢键数目会减少。

2.2.2   相对氢键能

抑制剂分子与水分子形成氢键的能量，与氢键数目和氢键的键长有关。氢键势能通常采用下式计算：

式中：U_hb为氢键势能，kJ/mol；$r$为受体原子与氢原子间的实际距离，即图4中的r_AH，10^–1 nm；D₀为氢键强度，kJ/mol；r₀为受体原子与氢原子间的平衡距离，10^–1 nm。

抑制剂分子与水分子所形成氢键的键长是不固定的，因此分析中采用了统计数据的平均值。抑制剂分子和水分子形成氢键的相对氢键能可表示为：

式中：H_en为形成氢键的相对氢键能；$k$为修正系数，与氢键强度有关，取1.0；$n$为形成氢键数量。

${r_0}$取平均值，利用式（3）计算了动力学抑制剂质量分数分别约为3.00%、8.50%和10.50%时，不同动力学抑制剂与水分子形成氢键的相对氢键能，结果如图6所示。

图  6  不同质量分数下各动力学抑制剂分子与水分子形成氢键的相对氢键能

Figure  6.  Relative hydrogen bond energies of hydrogen bonds formed by each kinetic inhibitor molecule with different mass fractions and water molecules

下载: 全尺寸图片 幻灯片

由图6可知，PVP，PVCap，PMC和YZ等4种动力学抑制剂在溶液中质量分数大体相同的情况下，抑制剂与水分子所形成氢键的相对氢键能随温度变化的趋势是不一致的。抑制剂质量分数为3.00%时，PVP分子和PMC分子与水分子形成氢键的相对氢键能随着温度升高而变大；YZ分子和PVCap分子与水分子形成氢键的相对氢键能随着温度升高而减小；抑制剂质量分数为8.50%时，PVP分子和PVCap分子与水分子形成氢键的相对氢键能随着温度升高缓慢地减小，而PMC分子和YZ分子与水分子形成氢键的相对氢键能随着温度升高迅速减小；抑制剂质量分数为10.50%时，PVP分子、PVCap分子和PMC分子与水分子形成氢键的相对氢键能随着温度升高均呈减小的趋势。抑制剂质量分数一定时，按其与水分子形成氢键的相对氢键能大小对动力学抑制剂进行排序，依次为PMC，YZ，PVP和PVCap。

3.   结　论

1）具有环状结构的PVP和PVCap可以有效降低水分子的扩散系数，从而减缓天然气水合物晶胞因温度升高导致的分解速率上升的趋势。

2）具有链状结构的PMC分子与水分子可以更好地结合形成氢键，且氢键的相对键能较大。

3）在所研究的PMC、YZ、PVP和PVCap等4种抑制剂中，PVCap比PVP有更好的温度适应性，即PVCap与水分子的相互作用受温度的影响更小，同时具有环状结构和链状结构的YZ其温度适应性最好。