TBAB半笼型甲烷水合物稳定性分子动力学模拟

2021-03-17梁海峰吉梦霞

天然气化工—C1化学与化工 2021年1期

张强，梁海峰，吉梦霞，杨琨

（太原理工大学化学化工学院，山西太原 030024）

随着常规化石能源的逐渐枯竭，人类对新能源的需求量越来越大，天然气水合物作为一种非常规天然气，资源总量非常庞大，估计为3.1 × 1015～7.6 ×1018m3，因此其开发利用得到了广泛的关注。天然气水合物是在低温、高压条件下烃类气体分子和水形成的类冰状晶体化合物，主体分子（水分子）通过氢键构成笼型骨架，客体分子（烃类分子）填充在晶穴中，主客体分子间的作用力为范德华力。

目前发现的水合物晶体包括Ⅰ型、Ⅱ型和H型，自然界中存在的水合物以Ⅰ型和Ⅱ型为主，但大约只有33%的晶穴被客体分子占据[1]。 Saim等[2]指出甲烷水合物在纯水体系中的客体分子晶穴占有率低于22.2%；Saruprias等[3]模拟了CO2水合物的分解过程，指出温度和客体分子晶穴占有率是影响水合物稳定性的重要因素；Sugahara等[4]指出了晶穴占有率对水合物稳定性影响的规律；Yu等[5]研究表明水合物笼中CH4浓度越高，水合物由亚稳态结构向稳定结构的转变速度越快；芦文浩等[6]研究了不同促进剂下甲烷水合物的稳定性，表明温度较低、客体分子与晶穴的直径比接近1时，水合物稳定性较好；耿春宇等[7]采用分子动力学方法模拟Ⅰ型甲烷水合物，指出由于客体分子占据水合物晶穴，晶穴间的相互作用能降低，水分子间的氢键作用加强，水合物的稳定性提高，并且在晶穴占有率高于87.5%的情况下，水合物表现出较好的稳定性。由此表明，较低的温度和较高的晶穴占有率能提高水合物的稳定性。

Fowler等[8]发现较大的客体分子破坏水合物笼型结构后，占据多个半笼型晶穴，形成特殊的水合物，称为半笼型水合物。四丁基溴化铵（TBAB）水合物是最常见的一种半笼型水合物，由于其具有良好的可降解性和低毒性，能有效加快反应速度和降低反应条件，且不容易挥发，因此受到广泛关注和研究。 Shimada等[9]发现TBAB水合物存在A型（TBAB·26H2O）和B型（TBAB·38H2O）两种结构，512晶穴能捕获气体小分子（CH4、N2、CO2等），形成气体水合物。Arjmandi等[10]和Mohammadi等[11]测量了H2/N2/CH4/天然气+ TBAB水溶液的相平衡数据。 Zhen等[12]研究表明，由于TBAB半笼型水合物只有512晶穴能捕获客体分子，因此TBAB半笼型水合物的选择性优于Ⅰ型、Ⅱ型和H型水合物；Aladko等[13]指出，TBAB水合物中最稳定的水合数为38H2O。目前，TBAB半笼型水合物的相平衡实验已得到充分研究，但由于晶体结构的复杂性，其微观结构的动力学内容尚未有研究，且水合物客体分子的晶穴占有率在实验过程中难以人为控制，因此本文基于分子动力学模拟，从微观角度分析了B型TBAB水合物的稳定性机理。

1 模型构建和模拟方法

1.1 模型构建

TBAB·38H2O水合物单晶胞属于正交晶系，空间群为pmma，水分子中氧原子和TBAB中原子的初始坐标根据XRD结果确定[14]，水分子和溴离子构成水合物笼，氢原子无序排列，且符合Bernal-Fowler规则[15]。模拟体系由2 × 2 × 2个TBAB半笼型水合物单晶胞构成，在x、y、z坐标方向上采用周期性边界条件，晶胞尺寸为4.2120 nm × 2.5286 nm × 2.4036 nm，体系包含608 个水分子，16 个四丁基氨阳离子（TBA+）和16个溴离子（Br-），48个512晶穴、32个51262晶穴和32个51263晶穴，其中TBA+填充在由51262晶穴和51263晶穴形成的半笼晶穴中，512晶穴由CH4分子填充。本文通过分子动力学模拟研究不同晶穴占有率（θ=100%、75%、50%、25%）在不同温度下对TBAB半笼型水合物稳定性的影响，模型初始结构如图1所示。

1.2 模拟方法

运用Materials Studio 8.0软件，采用COMPASSⅡ力场（Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies）在正则系综（NVT）下进行模拟，使用Ewald加和法描述长程静电作用力和范德瓦尔斯作用力，SPC（Simple Point Charge）模型描述水分子间相互作用，其中水分子视为刚性分子，键长和键角的值固定，Lennard-Jones势能函数描述水分子和其它分子间非键相互作用。

在进行分子动力学模拟之前，先对模型采用最速下降法和共轭梯度法进行几何优化，使结构处于能量最小状态。然后将溴离子和氧原子的坐标固定，使用Nose-Hoover方法控制体系温度分别在220 K、240 K、260 K、280 K和300 K下进行100 ps（1 ps=10-12s）的NVT驰豫，目的是完成结构松弛，减小结构误差，时间步长为1 fs（1 fs = 10-15s），每1000步输出一帧结构图；然后解除溴离子和氧原子的固定，在相同条件下，选择NVT系综进行动力学模拟，模拟时间为400 ps。图2为分子动力学计算过程中体系能量变化图，可以发现，在模拟的前100 ps内，体系的能量稳定在一个固定的值附近，说明此时体系已经达到了稳定状态，后300 ps可以分析体系的相关性质。对于每个体系，分子动力学模拟过程均进行三次，保证模拟的准确性。

2 结果与讨论

2.1 构象分析

模型的初始结构是理想的晶体结构，笼型结构框架具有清晰的氢键网络结构和良好的对称性（图

1），在分析笼型结构的稳定性时，可以分析氢键的规整程度和氧原子的对称程度来判断水合物笼型骨架的变形程度，进而直观的判断水合物的稳定性情况。

从图3可以看出在220 K时，体系A和B的水合物氢键结构变化很小，客体分子仍然位于晶穴的中心位置，水合物笼型完整，体系C和D的水合物笼型结构发生扭曲，部分客体分子偏离中心位置，此时水合物的稳定状态相比A和B体系下降，稳定性情况为A ＞B ＞C ＞D，由此得出晶穴占有率越高，水合物越稳定。图4中在相同占有率情况下，随着温度的升高，水合物笼型结构逐渐被破坏，水分子和TBAB逐渐偏离原始位置，且温度越高，偏离的位移越大。240 K时CH4分子开始从笼中逸出，说明水合物开始分解；260 K时体系的笼型结构开始坍塌，水合物两侧的CH4分子开始逸出，中心位置的晶穴仍能保持笼型结构，表明水合物的分解是由外侧向内侧逐渐演变的；到300 K时，体系的水合物笼型完全坍塌，CH4分子发生聚集，水合物完全分解，并出现分相现象。分析图4可知，在相同占有率的情况下，随着温度升高，水合物稳定性下降。

图1初始构象中水合物完整的笼型个数为48。分析不同温度下ABCD四个体系的最终构象图，可以得出表1，表中的数据为各构象图中完整的水合物笼的个数。纵向对比表中的数据，可以发现温度一定时，晶穴占有率越高，完整的笼型个数越多，水合物越稳定；同样的，当晶穴占有率相同时，随着温度的升高，完整笼型的个数减少，水合物稳定性下降。综合表1的分析也会得出，晶穴占有率越高，温度越低，水合物越稳定。

表1 各体系在不同温度下水合物稳定性情况（笼型个数）

2.2 径向分布函数

径向分布函数（Radial Distribution Function，RDF）是分析水合物结构特征最常用的评价参数，是晶体结构无序化的一个衡量指标，gαβ(r)表示在距离α粒子质心r处出现β粒子的概率，相当于体系的区域密度与平均密度的比值。当r值较小时，g(r)若出现几个极大值，则在这些r值的附近，出现β粒子的几率较大，此时结构处于有序状态；当r值较大时，g(r)值接近1，说明此时是无规则状态。若水合物处于稳定状态，则水分子中氧原子的排列具有周期性，RDF曲线有多个峰值，峰高随着r的增大逐渐降低；若水合物处于失稳状态，第二特征峰及其以后的特征峰将降低或消失。对于采用周期性边界条件的体系，RDF的定义如公式（1）所示：式中：V为模拟体系的体积，10-30m3；Nα和Nβ分别为α粒子和β粒子的个数；niβ(r)表示以第i个α粒子为中心，在r～r+Δr（单位为10-10m）范围内β粒子的数目。

图5为不同体系在不同温度下的水分子中氧原子的径向分布函数（RDF）。从图中可以看出，RDF曲线在r= 2.73 × 10-10m处出现第一特征峰，表明最相近的水分子之间氧原子的距离为2.73 × 10-10m，第二特征峰出现在r= 4.52 × 10-10m，这与文献[16]一致，验证了模型的正确性。随着温度的升高，峰值降低，峰谷升高，说明水分子无序性增大，笼型框架发生扭曲，水合物稳定性下降。在A和B体系中，温度升高到280 K时，第二特征峰消失，之后r值在1上下波动，与液态水的特征峰相似[17]，得出水合物开始分解。

对比在相同温度下不同体系的径向分布函数可以发现，CH4含量越大，RDF的特征峰越明显，峰值相对较高，峰谷较低，表明水分子有序性较强，水合物稳定性较好。 A和B体系差别不明显，说明晶穴占有率在75%以上，水合物处于相对稳定状态，此时温度对水合物的稳定性占主导因素。通过不同角度对比可以表明，温度越低，晶穴占有率越高，水分子的有序性越强，水合物越稳定。

2.3 均方位移

在分子动力学运动过程中，体系中的原子时刻在自由移动，每个瞬间各原子都会偏离原始位置到达某一位置。原子位移平方的平均值就称为均方位移（Mean Squared Displacement，MSD），如公式（2）所示。

式中：N为体系中的粒子总数；Ri(t)、Ri(0)分别表示i粒子在0、t时刻的位置。

均方位移反映了t时刻体系中原子的空间位置与初始位置的偏离程度。稳定的水合物，其分子处于相对固定的位置，MSD的值会在一个较低的值附近波动。若水合物发生分解，其分子位置坐标相对自由，MSD值会随着时间的增加而逐渐递增。

图6为各体系在不同温度下水分子中氧原子的均方位移图，A和B体系中，温度在260 K以下时，MSD值很小，表明水分子的位置相对固定，水合物笼型结构较完整，水合物处于稳定状态。随着温度进一步升高，MSD迅速增大，此时水合物开始分解，水分子处于无序状态，偏离初始位置较远，水合物笼型结构坍塌，且温度越高斜率越大，由此得出温度越高，水合物分解速率越快。

相对于A、B体系，C、D体系中只有220 K时MSD较小，说明低温有利于水合物稳定存在，且晶穴占有率越低，MSD越大，笼型结构变形程度较大，水合物极易分解，原因是晶穴占有率越低，水合物中的空笼越多，由于缺少客体分子的支撑，空笼的稳定性会下降，水合物的稳定性降低。对比四个图发现温度在260 K以下时，C、D体系水合物稳定性远远低于A、B体系，此时晶穴占有率对水合物稳定性的影响大于温度对水合物稳定性的影响。

2.4 氢键的数目与键长

氢键是一些与电负性较大的原子成键的氢原子和附近电负性较大或者带孤对电子的原子形成的较强的非键作用，其作用强度介于成键作用和非键作用之间。体系中的氢键包括水分子之间的氢键和溴离子与水分子之间形成的氢键。当水合物处于稳定状态时，氢键数目与长度相对于初始结构的氢键变化较小，若水合物发生分解，则氢键的数目与长度变化较大。本文所选的COMPASSⅡ力场对氢键的统计是基于对范德华力和静电相互作用的综合分析而实现的。

与RDF和MSD不同，氢键的变化是描述水合物笼型结构随温度和晶穴占有率变化而变化的直观体现。水合物初始结构中，氢键的数目为1215，平均键长为1.82 × 10-10m。在分子动力学过程中，随着水合物的分解，部分氢键断裂，部分氢键被拉长，导致水合物笼型扭曲。图7为分子动力学之后水合物中氢键的数目和平均键长的变化，从图中可以看出，随着温度升高，氢键数目不断减少，键长逐渐变长，表明水合物笼出现不同程度的扭曲变形；在同一温度下，晶穴占有率越低，氢键数目越少，键长越长，水合物越不稳定。通过分析氢键的变化也可直观地看出低温和高占有率有利于水合物的稳定存在。