分子模型SPC/Fw-EPM2对CO2水合物动力学及热力学特性表征

苏保照 焦丽君,2 李裕斌 李 倩

（1.青岛职业技术学院海尔学院 青岛 266555；2.中国石油大学新能源学院 青岛 266555）

0 引言

CO2水合物由水和二氧化碳分子组成，水分子之间通过氢键相连形成笼形结构，二氧化碳分子位于笼子中心，依靠范德华力相互作用[1]。CO2水合物能够在0-10℃温度条件下生成，与空调运行温度相适宜，而且相变潜热较大，约为500kJ/kg，远大于冰的分解焓，是一种储能性能良好的相变材料，未来极有可能取代冰蓄冷，解决空调储能需要双工况机组的技术难题[2,3]。而且与N2、H2和CH4等水合物相比，相平衡条件更加温和，用于烟气中二氧化碳捕集，能够解决传统碳捕集技术费用过高的问题[4-6]，为碳中和、碳达峰目标的实现提供技术支撑。

然而，由于CO2水合物合成过程存在诱导期长、过冷度大和生长缓慢等问题，限制了相关技术应用。为此，研究人员通过搅拌、鼓泡等动态强化技术和使用表面活性剂、纳米粒子等静态强化手段促进水合物合成，发现实验条件对动力学促进效果影响较大，已有实验研究中关于CO2水合物合成的动力学性质差别较大，而且作用机理尚不明确[7]。对于水合物分解动力学过程，研究人员发现传热方式、加热温度、降压速度、添加剂种类和浓度等对水合物分解速率影响较大，但是分解过程动力学方程大多只能较好的拟合特定实验，不具有普遍适用性[8-12]。由于实验条件的限制，只能观察实验现象和宏观动力学规律，若要洞悉水合物动力学行为背后的非平衡态传热传质规律，还需要从分子水平进一步研究。

分子动力学模拟作为一种微观研究手段，能够从分子尺度表征微观特性，但是对分子模型的依赖性较强[13]。水分子模型主要有刚性分子模型SPC/E、TIP4P/Ice 和TIP4P/2005，柔性分子模型SPC/Fw，二氧化碳分子模型包括TraPPE 和EPM2[14]。其中，柔性分子模型SPC/Fw 在键长伸缩和键角弯曲两个方面具有自由度，相比于刚性分子模型，对温度的描述更加准确，分子模型作用下的热力学性质或许也更加准确。为此，使用分子模型SPC/Fw 对CO2水合物的微观特性进行表征十分有必要。

1 模拟方法

1.1 模拟体系

为模拟CO2水合物热力学及动力学特性，运用Material Studio 软件建立水-水合物-二氧化碳三相体系，如图1 所示。左侧区域包括368 个水分子，采用棒状模型，深色表示氧原子，浅色表示氢原子，氢键由虚线表示。右侧区域包括192 个二氧化碳分子，采用球状模型，深色表示氧原子，浅色表示碳原子。中间区域由2×2×6 超晶胞二氧化碳水合物组成，包括192 个二氧化碳分子和1104 个水分子，其中氧原子和氢原子按照Takeuchi 等报道的笛卡尔坐标放置，二氧化碳分子随机放置在笼子中心[15]。模拟体系沿x、y 和z轴的尺寸为24.06Å×24.06Å×128.72Å。分子模型选用SPC/Fw和EPM2，具体参数见表1 和表2[16,17]。分子间作用力包括范德瓦尔力、静电力和键长伸缩、键角弯曲作用力，如公式（1）所示，不同原子之间的势能参数采用Lorentz-Berthelot 规则计算，如公式（2）所示[18]。

图1 模拟体系初始时刻构型Fig.1 Initial configuration of the simulated system

表1 势能参数Table 1 Potential energy parameters

表2 键参数Table 2 Bond Parameters

式中，q为电荷，ε0为真空介电常数，r为原子之间的距离，ε和σ为势能参数，k为弹力常数，r和θ分别为键长和键角。

1.2 模拟过程

分子动力学模拟运用开源软件LAMMPS 进行，其中data 文件由模拟体系转化成坐标数据得到，in 文件设置参数为：采用周期性边界条件；分子间相互作用力采用长程静电力，截断半径为12Å，计算方法采用pppm（particle-particle particle-mesh），计算精度取10-6；时间步长为1fs，使用velocity-verlet 积分算法求解牛顿运动方程；温度和压力控制使用Nose-Hoover恒温器和恒压器处理，阻尼系数分别设为0.1ps 和1ps。模拟步骤为：首先使用共轭梯度算法对体系进行能量最小化，然后将体系置于NPT 系综下弛豫100ps，使体系达到目标压力和温度，最后将体系置于NVE 系综下运行约5ns。

2 结果与讨论

2.1 微观结构特性

径向分布函数表示局部密度与体系密度的比值，可以表征体系的微观结构，定义如公式（3）所示。

式中，N为分子总数，r为距离“参考分子”中心的距离，ΔN为介于r到r+δr之间分子数，ρ为系统密度，δr为设置的距离差，T为计算的总时间。水中氧原子间径向分布函数用go-o(r)表示，二氧化碳中碳原子间径向分布函数用gc-c(r)表示，结果如图2 所示。可以看出，对于水中氧原子间径向分布函数，当r=2.73Å时，go-o(r)最高，即处于第一个峰值，与水合物中氢键形成的氧原子距离一致，当r=4.47Å和6.45Å时，处于第二和第三个峰值，与水合物中四面体网络中氢键结构以及六元环中氧原子间距离相一致。而且与Kondori 等[13]报道的结果基本一致（三个峰值分别为2.75、4.5 和6.53Å）。对于二氧化碳中碳原子间径向分布函数，当r=6.63Å时，gc-c(r)最高，表示笼子中二氧化碳分子之间的距离，由于分子振动或转动，峰值在一定范围内波动（6.45～6.87Å），与文献报道一致[19]。这表明使用分子模型SPC/Fw-EPM2 能够准确描述CO2水合物的微观结构特性。不止是径向分布函数，由LAMMPS 软件统计输出沿z 轴的分子密度分布也能证明模型的可靠性，如图3 所示，二氧化碳和水分子均有11 对高低峰，与图1 所示体系构型一致，液态水区和二氧化碳集中区域也呈现出质量密度不变的特点。

图2 初始时刻不同原子间径向分布函数Fig.2 The radial distribution function between different atoms at the initial moment

图3 初始时刻不同分子沿z 轴质量密度分布Fig.3 Mass density distribution of different molecules along the z-axis at the initial moment

2.2 分解动力学特性

本征动力学理论认为水合物在分解过程中存在活化能垒，可以作为定量描述水合物动力学特性的基本参数。计算方法如公式（4）和（5）所示[20]。

为拟合不同温度条件下CO2水合物分解速率，运用SPC/Fw-EPM2 分子模型分别在3.5MPa、300-320K 初始热力学条件下模拟，类水合物二氧化碳分子数目由Myshakin 等[22]报道的判别方法计算，数目随时间变化如图4 所示。可以看出，初始温度越高，水合物分解速率越快，这与普遍的实验结果一致。对不同温度条件下，类水合物二氧化碳分子数目随时间变化作线性拟合，得到分解速率，如表3 所示。将公式（4）和（5）合并，得到分解速率与温度的依赖性表达式及活化能计算方法，如公式（6）所示。

表3 温度与分解速率数据统计Table 3 Statistics of temperature and decomposition rate

图4 不同热力学条件下类水合物二氧化碳分子数目随时间变化Fig.4 Number of hydrate-like carbon dioxide molecules change with time under different thermodynamic conditions

由此将分解速率自然对数与温度倒数进行拟合，结果如图5 所示，拟合优度为0.9790，拟合斜率为-6093.82±514.94，求得活化能ΔE为50.66±4.28kJ/mol，与实验值在同一数量级（73kJ/mol）[23]，而且比前期使用SPC/E-EPM2 分子模型计算的活化能（18.48±0.45kJ/mol）更接近实验值[14]，一方面表明分子模型SPC/Fw-EPM2 在描述CO2水合物的动力学特性方面的适用性和准确性，另一方面，也体现了模拟体系和研究方法的正确性。

图5 分解速率对初始温度的依赖性Fig.5 Dependence of decomposition rate on initial temperature

2.3 热力学性质

尽管水合物的热力学特性已经获得了相对比较精确的研究，但是在模拟过程中由于热力学条件的改变，相变动力学特性差别比较大，这一方面可能与分解界面的性质有关，另一方面也与分子模型的稳定性有关，即不同分子模型表征的水合物相平衡条件差别很大。Costandy 等使用四力点模型水分子模型TIP4P/Ice 和TIP4P/2005 及二氧化碳分子模型TraPPE，对CO2水合物相平衡特性进行了模拟和预测，尽管取得了较好的实验一致性，但是由于TIP4P/Ice 和TIP4P/2005 均属于刚性分子模型，在需要温度梯度和显热热流描述的模拟研究中并不适用。水合物的相平衡特性可以根据体系势能随时间变化的行为确定，对于图1 所示三相模拟体系，若保持压力不变，当温度高于相平衡温度时，水合物分解，分子间随着距离的增大相互作用力减小，表现为势能上升，反之，新相水合物合成放热，势能下降。但是由于合成过程成核所需诱导期长，受到计算资源限制，仅模拟5ns 时间，故仅针对体系的热力学稳定性作相关说明。以3.5MPa、275-290K初始条件为例，模拟体系的势能和温度变化如图6所示。可以看出，290K 初始温度下，体系势能显著升高、温度显著下降，说明水合物分解。275K初始温度下，体系势能和温度均不变，说明水合物处于稳定区。278 和280K 初始条件下，体系势能在分解前期缓慢上升，随着温度的下降，当降到平衡温度附近时，没有足够的温差为分解过程提供动力，而且随着二氧化碳气体的释放，体系压力增大，水合物分解停止。从图中看，平衡温度约在275K左右，比相似压力下测得的实验值略低（3.395-3.667MPa、280.34-280.79K）[24]。

图6 不同初始条件下（a）模拟体系势能和（b）温度随时间变化Fig.6 (a)Potential energy and(b)temperature of the simulated system change with time under different initial conditions

Baghel 等[25]提出平衡温度可以使用公式（7）计算[24]。

式中，T为瞬时温度，Teq为平衡温度，T0为初始温度，k为拟合系数。以310K 初始温度条件为例求解平衡温度，拟合结果如图7 所示，得到平衡温度为279.85±0.33K，与实验结果比较接近。正是由于柔性分子模型SPC/Fw 在键长伸缩和键角弯曲两个方面具有自由度，相比于刚性分子模型，对温度的描述更加准确，分子模型作用下的热力学性质也更加准确。这一点在Li 等[26]对甲烷水合物的研究中也得到了证实。对于310K 初始温度，水合物在4000ps 时已分解结束，体系压力增大，因此平衡温度也比275K 略大。

图7 310K 初始温度条件下模拟体系平衡温度Fig.7 Equilibrium temperature of simulated system at the initial temperature of 310K

使用Clausius-Clapeyron 方程求解分解热，如公式（8）所示，相平衡温度和压力关系取自文献[24]，由公式（7）拟合后得到斜率为-10323。

式中，ΔHd为分解热；R为通用气体常数，大小为8.314J·mol-1·K-1；p为压力，MPa；T为温度，K；Z为压缩因子，利用PR 方程计算[27]，大小为0.6954，分解热为59.68kJ·mol-1，在文献报道的数值范围内（54.5-60.1kJ·mol-1）[28]。

此外，建立纯水合物结构体系，运用平衡分子动力学模拟方法，借助Green-Kubo 公式，编程计算了水合物热导率。计算公式如下：

式中，kB为波尔兹曼常数，大小为1.38×10-23J/K；V为模拟区域的体积，m3；T为模拟区域的温度，K；为热流自相关函数（HCACF）。270K、3.5MPa 条件下，热导率随时间变化如图8 所示，后期热导率趋于稳定，平均热导率为0.84W·m-1·K-1，在Jiang 等[29]报道的二氧化碳水合物热导率范围内（0.75-0.95）W·m-1·K-1，表明分子模型SPC/Fw-EPM2 能够较好的描述水合物的热力学性质。

图8 270K、3.5MPa 条件下热导率随时间变化Fig.8 Thermal conductivity changes with time at 270K,3.5Mpa

3 结论

考虑到柔性分子模型在温度描述中的精确性，运用分子模型SPC/Fw-EPM2，在不同热力学条件下，通过分子动力学模拟对CO2水合物相变过程动力学行为和热力学特性进行了微观表征，从微观结构特性、分解活化能和热力学稳定性三个方面表述了分子模型的适用性，为基于温度梯度和显热热流统计的涨落耗散分析提供了理论基础，为水合物动力学行为的微观机理探索提供了理论支撑。具体结论包括：

（1）分子模型SPC/Fw-EPM2 作用下，水中氧原子间径向分布函数的第一、第二和第三峰值分别为2.73、4.47 和6.45Å，二氧化碳中碳原子间径向分布函数的峰值在6.45-6.87Å之间波动。分子模型SPC/Fw-EPM2 能够准确描述CO2水合物的微观结构。

（2）3.5MPa、300-320K 初始热力学条件下，CO2水合物分解速率对初始温度的依赖性较强，活化能ΔE为50.66±4.28kJ/mol，分子模型SPC/Fw-EPM2 能够准确描述CO2水合物的动力学特性。

（3）3.5MPa、270K 条件下CO2水合物处于热力学稳定区域，3.5MPa、290K 条件下CO2水合物分解。310K 初始温度条件下，模拟体系平衡温度为279.85±0.33K，分解焓为59.68kJ·mol-1。3.5MPa、270K 条件下CO2水合物的平均热导率为0.84W·m-1·K-1，均与实验结果相接近，表明分子模型SPC/Fw-EPM2 能够准确描述CO2水合物的热力学性质。