廖清云，史博会，宋尚飞，王君傲，姚海元，段 旭，宫 敬，闻 峰

（1. 中国石油大学（北京） 天然气水合物国家重点实验室，油气管道输送安全国家工程实验室，石油工程教育部重点实验室，城市油气输配技术北京市重点实验室，北京 102249；2. 中海油研究总院有限责任公司，北京 100028；3. 国家石油天然气管网有限公司，北京 100007）

随着深海油气资源开发的快速发展，气体水合物的开发利用及其在流动保障中的防治研究逐渐深入，其中，气体水合物生长及分解过程的微观机理成为了研究热点。考虑到常规实验手段在时间尺度和空间尺度上的限制，国内外已有多位学者利用分子动力学模拟进行了微观尺度的研究，包括水合物在纯水体系中的生长机理[1]、添加剂或杂质（如无机盐/水合物动力学抑制剂/天然果胶/蜡分子等）的存在对水合物生长机理的影响[2-5]、水合物在逐步升温过程中的分解机理[6]、水合物在无机盐中的分解机理[7]及不同温度和压力条件下抑制剂存在时水合物的分解机理[8]等。

在水合物生长及分解过程的分子模拟研究中，均需在初始模型中添加根据实验数据[9]搭建的水合物晶胞结构作为水合物生长种子或水合物分解前体结构。因水合物晶胞结构中氢原子位置具有不确定性，在分子模拟中需根据Bernal-Fowler规则、体系构型的势能最小原则或氢键水分子朝向成比例原则等，进一步确定氢原子坐标。在构建模拟模型后，需开展能量最小化以降低体系能量并寻找体系势能极小值以消除体系内不合理接触，随后通过弛豫过程，即一段有限时长的分子动力学模拟，以进一步调整能量最小化后的模型，使系统从初始配置逐渐演变达到平衡状态，才能为分子模拟计算提供更为合理的初始构型，该过程对获得合理的分子模拟结果至关重要。

弛豫过程以分析体系的热力学参数（如温度、能量和密度等）及结构参数（如序参数、均方位移和径向分布函数等）为主[10]。然而，不同研究者选择的弛豫系综及弛豫时间各有不同：Tung等[1]以具有固定温度的正则系综（NVT）（T= 200 K）弛豫20 ps平衡气-液-固体系；Li等[3]以NVT系综（T= 276 K）弛豫400 ps平衡液-固体系；Xu等[4]以NVT系综（T= 260 K）弛豫30 ps，及具有固定温度及压力的系综（NPT）（T= 260 K，p= 15 MPa）弛豫30 ps平衡液-固体系；Choudhary等[6]以NPT系综（p= 100 MPa）弛豫5 ns平衡水合物晶胞结构；Xu等[7]以NVT系综在模拟温度弛豫50 ps平衡液-固体系；Kondori等[8]以NVT系综（T= 273 K）弛豫50 ps平衡水合物晶胞结构；张庆东等[11]首先以NVT系综（T= 263.15 K、T= 273.15 K、T= 283.15 K）固定氧原子位置弛豫50 ps，随后解除氧原子约束在相同条件下弛豫1 ns平衡水合物晶胞结构。总之，NVT系综在弛豫中被广泛应用，且弛豫时间通常设置在皮秒量级，但关于如何确定合理的弛豫方案的研究较少。

本文以甲烷（CH4）水合物晶胞笼型结构的演变规律为分析参数，探讨不同弛豫方案对CH4水合物晶胞构型的影响，包括：弛豫系综、弛豫时间、弛豫压力与温度等，目的是确定更为合理的弛豫方案以开展CH4水合物的分子动力学模拟，为水合物分子动力学模拟研究提供理论基础。

1 模拟模型和方法

1.1 CH4水合物晶胞模型

图1 3 × 3 × 3 sI型CH4水合物晶胞示意Fig. 1 Schematic diagram of 3 × 3 × 3 sI methane hydrate unit cell

1.2 模拟设定

基于相关文献研究[1,3,5,16]，对于分子间相互作用，确定选用TIP4P-Ew力场[17]描述H2O分子，OPLS-AA力场[18]描述CH4分子。各原子的电荷及L-J势函数参数如表1所示，不同原子间的L-J势函数参数值通过Lorentz-Berthelot混合规则计算。将H2O分子视为刚体，采用SHAKE算法[19]约束H2O分子的H-O键长与H-O-H键角。范德华力采用L-J势函数计算，截断半径为1 nm。短程库伦作用截断半径为0.9 nm，长程库伦相互作用通过pppm算法计算[20]，计算精度为10-4。模拟时间步长为1 fs，每0.5 ps输出一次热力学参数与原子位置信息。利用Nose-Hoover恒温器及恒压器实现体系的温度压力控制，模拟体系在x、y和z方向上均采用了周期性边界条件。

表1 L-J势函数参数值与原子电荷值Table 1 L-J potential function parameter value and atomic charge value

1.3 CH4水合物晶胞构型稳定性判定方法

本文以一种改进的笼型结构识别方法判定CH4水合物晶胞构型稳定性，其具体流程如下：

（1）筛选成笼氧原子。当相邻氧原子间的间距小于0.35 nm，且相邻的两个氧原子及与其相连的一个氢原子间形成的角度小于30°，则认为这两个氧原子之间形成氢键[21]，当一个氧原子上形成的氢键数大于等于3时则认为其为成笼氧原子。

（2）基于氧原子筛选结果，进行五元环的识别。将成笼氧原子视为无向图的顶点，根据各氧原子之间的连接关系来判断是否形成五元环。

（3）512笼结构如图2所示，绿色实线为512半笼的中心五元环，红色实线为半笼中与中心五元环相连的五元环，蓝色实线为完整的512半笼结构。考虑到其为对称结构，可采取半笼拼接的方法进行识别[22]。根据五元环的识别结果，判断各五元环之间的连接关系，进而确定以各五元环为中心的512半笼，再将得到的512半笼进行拼接，最终得到512全笼。

图2 512笼结构示意Fig. 2 Schematic diagram of 512 cage structure

水合物晶体构型内氢原子的位置会影响氢键判定中的角度计算结果，并影响所识别的环笼结构数目。因此，以该改进的笼型识别方法分析弛豫过程中体系笼型结构数目，若所识别的笼型结构数目越接近理论值且波动越小，则说明该弛豫方案的设定越优。对于3 × 3 × 3 sI型CH4水合物晶胞结构而言，按理论分析应有35个512笼，其连接关系如图3所示，其中红色小球为氧原子，蓝色实线为氢键。

本文的成本系数(Ck)是指该方案的生命周期成本(工程造价)与所有方案成本总和之比。各方案的价值系数Vk = Fk / Ck，计算结果见表3。可见在4种防腐方案中，D方案的价值系数最高，即采用功能性玻璃钢复合材料作为火电厂烟囱防腐工程材料是最优解决方案。

图3 3 × 3 × 3 sⅠ型CH4水合物晶胞中512笼的正交示意(a)及透视示意 (b)Fig. 3 Orthogonal schematic diagram (a) and perspective schematic diagram (b) of 512 cages in 3 × 3 × 3 sⅠmethane hydrate unit cell

2 结果与讨论

2.1 弛豫系综及弛豫时间对CH4水合物晶胞构型的影响

为研究弛豫系综及弛豫时间对CH4水合物晶胞构型的影响，设定弛豫方案1-1及弛豫方案1-2。弛豫时间均为500 ps；弛豫方案1-1设定为NVT系综，弛豫方案1-2设定为NPT系综；以水合物-液态水-甲烷气三相平衡温压条件为基础，10 MPa时的水合物-液态水-甲烷气三相平衡温度约为281 K[1]，故暂设定弛豫温度为250 K，弛豫压力为10 MPa。为避免瞬时结构波动对结果分析造成影响，将512笼识别数量在3 ps的时间窗口内进行平均后再进行分析。

方案1-1（NVT，T= 250 K，下同）弛豫500 ps的512笼识别结果如图4(a)所示。在500 ps弛豫过程中，512笼数量变化可分为3个阶段，如表2所示。在58 ps内快速达到近平衡状态，在58 ps后随着模拟时间的延长逐渐趋于更稳定的平衡状态。此外，体系温度在± 5%范围内波动，如图4(b)所示，表明体系温度已达到平衡状态。

表2 不同弛豫方案下512笼数量变化阶段Table 2 Change stages of 512 cage number under different relaxation schemes

图4 方案1-1弛豫500ps的512笼识别结果(a)及温度变化(b)Fig. 4 512 cage recognition results (a) and temperature changes (b) of scheme 1-1 relaxation 500 ps

方案1-2（NPT，T= 250 K，p= 10 MPa，下同）弛豫500 ps的512笼识别结果如图5(a)所示。在500 ps弛豫过程中，512笼数量在弛豫初期波动较大，其变化可分为4个阶段，如表2所示。弛豫在20 ps内快速达到近平衡状态，在20 ps后随着模拟时间的延长逐渐趋于512笼数量更高的近平衡状态，在364 ps后达到更稳定的平衡状态。此外，体系温度在± 5%范围内波动，如图5(b)所示，表明体系温度已达到平衡状态。

图5 方案1-2弛豫500 ps的512笼识别结果(a)及温度变化(b)Fig. 5 512 cage recognition results (a) and temperature changes (b) of scheme 1-2 relaxation 500 ps

由表2可知，500 ps弛豫结束后，两种弛豫方案下体系中识别得到的512笼数量难以稳定达到35个。这是由于H2O分子在模拟中不断运动，其原子存在越过模拟盒边界的情况，进而导致边界处H2O分子所形成的氢键数未达到成笼条件而被忽略。方案1-1和1-2弛豫后的体系结构如图6所示，绿色实线为氢键，模拟盒边缘的各方向应各有3个512笼，因本环笼识别代码中的氢键筛选过程暂未考虑周期性边界条件，致使方案1-1未识别出的512笼位于模拟盒的左下边界及右中边界，方案1-2未识别出的512笼位于模拟盒的右下边界，故512笼数量最终呈现出在接近理论值附近波动变化的现象。

图6 方案1-1(a)及方案1-2(b) 弛豫后的体系结构Fig. 6 512 cages of scheme 1-1 (a) and scheme 1-2 (b) after relaxation

对于弛豫系综，对比方案1-1和方案1-2得到的512笼识别结果，如图7所示。由图7(a)可知，在方案1-2下，512笼数量能更快达到较高的数值，但其在364 ps内的数值波动较大。经对比发现，在较短的模拟时长内，方案1-1得到的结果稳定性更好。因此，可以认为方案1-2（NPT系综）下弛豫对笼型结构的形成有一定的促进作用，但在方案1-1（NVT系综）下弛豫更利于结构的稳定。这一现象与Zhang等[23]的结论一致：在NPT系综下水合物模拟的成核速率最快，而在NVT系综下模拟得到的水合物结晶度更高。由图7(b)可知，方案1-1和1-2弛豫过程中体系温度变化相近，均在平衡值附近波动。

图7 方案1-1与方案1-2弛豫的512笼识别结果(a)及温度变化(b)对比Fig. 7 Comparison of 512 cage recognition results (a) and temperature change (b) between scheme 1-1 and 1-2

对于弛豫时间，由图4(a)与图5(a)可知，弛豫时间越长，越有利于获得稳定的CH4水合物晶胞构型，但弛豫时间的增加会显著延长计算时间。适宜的弛豫时间，不仅可以避免因弛豫时间较短不能使体系达到稳定平衡的缺点，还可以弥补弛豫时间较长而导致计算时长增加的缺陷。对于波动较大的弛豫方案1-2，由于512笼数量在模拟初期变化明显，需选取适当较长的弛豫时间以达到体系的平衡，考虑到其在100 ps后波动幅度有所减小，且此时两种弛豫方案内识别得到的512笼数量均可达到32个，故选用100 ps作为本文后续分析研究的弛豫时间。

基于方案1-1及方案1-2，为进一步研究弛豫温度和压力条件及系综联用对CH4水合物晶胞构型的影响，设定多种弛豫方案，如表3所示。设定弛豫时间为100 ps，分别设置不同温度和压力以分析弛豫温压条件对弛豫过程的影响，并研究NVT及NPT系综联合应用对弛豫过程的影响。

2.2 弛豫温度对CH4水合物晶胞构型的影响

表3中弛豫方案2-1至2-3，是在弛豫方案1-1基础上，调整弛豫时间为100 ps，改变弛豫温度设定的3个方案。方案2-1至2-3的弛豫温度分别设置为200 K、250 K和300 K。3种不同弛豫温度下的笼识别结果对比如图8(a)所示。对于最低弛豫温度200 K，在相同弛豫时间下，体系内识别得到的512笼数量最高，且其波动幅度最小。然而，对于最高弛豫温度300 K而言，体系内所识别的512笼数量波动幅度最大。弛豫温度为250 K时，512笼识别数量结果处于两者之间。总体而言，在弛豫时间为75~100 ps时，弛豫方案2-1至2-3所识别的512笼数量逐渐接近。

图8 方案2-1至2-3弛豫的512笼识别结果(a)及温度变化(b)Fig. 8 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 2-1 to 2-3

表3 弛豫方案Table 3 Relaxation schemes

因水合物的生成条件是低温高压，在低温条件下，水合物生成的驱动力更强，故较低的温度有利于水合物笼型结构的稳定，更利于笼型结构识别分析。但随着弛豫时间达到100 ps，不同弛豫温度所达到的水合物晶胞构型稳定效果逐渐接近。此模拟结果与张庆东等[11]得到的结论一致，即弛豫温度会影响CH4水合物的稳定性，且弛豫温度越低，其稳定性越高；但对于已稳定的水合物晶胞结构，弛豫温度的影响将不再明显。因此，图8(a)的模拟结果也再次佐证了本模拟体系内CH4水合物晶胞构型在弛豫100 ps时已基本达到稳定。方案2-1至2-3弛豫过程中的温度变化如图8(b)所示，可知各体系的温度基本在设定值附近波动。

越低的弛豫温度越有利于CH4水合物晶胞构型的稳定，而过高的弛豫温度会减弱CH4水合物晶胞构型的稳定性。对于TIP4P-Ew水分子模型，冰点约为245 K[24]，若弛豫体系中含有液态水，需设置高于冰点的模拟温度以避免冰的产生。对于本文所选参数，三相平衡温度在10 MPa时为281 K，弛豫温度需低于所选力场模型的相平衡温度以避免水合物结构失稳。对于本文研究的体系，250 K的弛豫温度较合理，为此，表3中的弛豫方案3-1、3-2、4-1、4-2和4-3均选取250 K作为弛豫温度，以开展后续研究。

2.3 弛豫压力对CH4水合物晶胞构型的影响

表3中弛豫方案3-1和3-2，是在弛豫方案1-2基础上，调整弛豫时间为100 ps，改变弛豫压力设定的两个方案。方案3-1与3-2的弛豫压力分别设置为10 MPa与30 MPa。方案3-1和方案3-2在100 ps内的512笼识别结果对比如图9(a)所示。数据表明，增大压力后512笼数量的变化更加剧烈。在100 ps模拟末，30 MPa弛豫压力下的512笼识别数量仍有显著波动。这是因为分子动力学模拟通过调整模拟盒尺寸并随之缩放各原子位置以实现压力的控制，在高压下体系的体积调整较大，故其结构的波动变化也越剧烈，具体表现为在高压下弛豫时，所识别的512笼数量稳定性差。方案3-1和方案3-2下体系温度的变化如图9(b)所示，数据表明体系温度基本在设定值附近波动。因此，对于NPT系综下弛豫压力的确定，一方面，要避免高压所引起的体系结构变化；另一方面，应满足后续正式模拟压力需求。为此，针对本文研究体系，为保证结构稳定，选取10 MPa作为NPT系综弛豫压力。

图9 方案3-1和3-2弛豫的512笼识别结果(a)及温度变化(b)Fig. 9 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 3-1 to 3-2

2.4 双系综弛豫方案对CH4水合物晶胞构型的影响

经前续模拟结果分析，在适宜的弛豫时间和弛豫温度下，NVT系综下的弛豫能得到稳定的结构，但水合物结构达到一定数量所需的时间更长；在适宜的弛豫温度和压力下，NPT系综下可更快得到一定数量的水合物结构，但水合物结构的波动变化较大；两个系综下的弛豫均可平衡体系温度。因此，考虑将两个系综相结合，借助NVT系综易于得到稳定结构的优势，加以NPT系综对笼型结构形成的促进作用，得到表3中的弛豫方案4-1，总弛豫时长为200 ps，首先在NVT下弛豫100 ps以稳定初始结构，再在NPT下弛豫100 ps以进一步达到接近理论结构的稳定状态。

方案4-1模拟过程中512笼识别数量的变化如图10(a)所示。其亦可大致分为3个阶段，如表4所示。弛豫在58 ps内快速达到近平衡状态，在58 ps后随着模拟时间的延长逐渐趋于更稳定的平衡状态。方案4-1弛豫过程中体系的温度变化如图10(b)所示，弛豫过程体系温度基本在设定值附近波动。因此，在弛豫方案4-1下，首先进行NVT弛豫再进行NPT弛豫可有效地平衡体系。

图10 方案4-1弛豫的512笼识别结果(a)及温度变化(b)Fig. 10 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 4-1

表4 方案4-1 512笼数量变化阶段Table 4 Change stages of 512 cage number of relaxation scheme 4-1

将方案4-1、4-2及4-3进行对比，NVT弛豫后再进行NPT弛豫，512笼数量在NPT弛豫初期会出现一定波动，但在167 ps即达到33个512笼，且数值基本稳定，这一阶段快于如图11(a)所示的单独NVT弛豫，且其在100~200 ps弛豫时间内波动范围明显小于如图11(b)所示的单独NPT弛豫。在200 ps的弛豫时间下，单独进行NPT弛豫识别得到的512笼平均值更大，单独进行NVT弛豫识别得到的512笼平均数量最少。以上结果进一步表明，NPT弛豫有利于笼型结构的形成，但其结构稳定性较差，而NVT弛豫可优化稳定水合物结构，两者结合后有利于水合物弛豫平衡。

图11 方案4-3 (a) 及方案4-2 (b) 与方案4-1弛豫的512笼识别结果对比Fig. 11 Comparison of 512 cage recognition results of relaxation between scheme 4-1 and scheme 4-3 (a) and scheme 4-2 (b)

此外，在NPT弛豫过程中，由于体系体积处于波动变化，故以分段密度变化考察体系平衡状态，对所选方案的弛豫结果进行验证。对于NPT弛豫过程，从体系基本稳定（512笼达到31个）后开始每隔10 ps取一次密度的平均值，计算每段平均值的标准差，当标准差小于0.001 g/cm3时[25]，认为体系达到平衡状态。对于单独NPT弛豫，从模拟20 ps后（表2）开始取密度平均值；对于NVT + NPT弛豫，从NVT弛豫结束后，即100 ps开始取密度平均值。密度计算结果如表5所示，单独NPT弛豫在200 ps内分段密度标准差仍大于0.001 g/cm3，到500 ps时分段密度达到平衡状态。而首先进行NVT弛豫，再进行NPT弛豫，可在200 ps时达到密度的平衡状态。因此，与水合物笼型结构数量的变化类似，从分段密度的角度衡量，将NVT系综及NPT系综结合进行弛豫，水合物结构更易达到平衡状态。

2.5 双系综弛豫后水合物分解温度的模拟分析

经不同弛豫方案的模拟比较，拟推荐采用双系综弛豫方案4-1，即首先进行100 ps的NVT弛豫，再进行100 ps的NPT弛豫，故对双系综弛豫后处于平衡态的水合物构型进行分解模拟以验证弛豫方案的适用性。采用空隙诱导法进行水合物分解模拟，Choudhary等[6]的模拟表明，水合物分解温度在空隙率为10%左右时处于平衡阶段，此时的分解温度可确定为该力场下的水合物分解温度。基于双系综弛豫后的水合物构型，随机移除H2O分子及CH4分子共计115个，获得空隙率约为8%的水合物结构，如图12所示，橙线为氢键，绿球为CH4分子中的碳原子，蓝线为随机移除H2O分子后破坏的氢键，蓝球为随机移除的CH4分子。

水合物分解模拟过程相关参数设置与弛豫过程一致，如1.2节所述。模拟压力为10 MPa，模拟温度由随机移除分子后的平衡温度185 K匀速升温至350 K，升温速率为0.165 K/ps，总模拟时长为10 ns。体系的温度及势能变化如图13所示，势能突变位置即为相变点，故该模型在10 MPa下的分解温度约为292 K。对比10 MPa下水合物实际分解温度286.2 K，模拟分解温度高约5.8 K。这一偏差是由于水合物分解温度模拟研究依赖于模拟体系大小，3 × 3 × 3 sI型CH4水合物晶胞模拟系统较小，而较小体系模拟所得到的水合物分解温度较高[6]。此外，Choudhary等[6]使用3 × 3 × 3 sI型CH4水合物晶胞进行分解模拟时，所得分解温度偏差范围为5~35 K，故本模拟结果在偏差范围内且偏差较小。对于水合物分解过程，某一水合物笼坍塌后，水合物开始分解，如图14(a)所示（6.5 ns时），在随后的1 ns内完全分解，由于分解后CH4分子过饱和，形成气液两相体系，如图14(b)所示（7.3 ns时），其中橙线为氢键及H2O分子，绿球为CH4分子中的碳原子。以上模拟结果表明，双系综弛豫方案所得的稳定水合物构型可用于后续产生相的模拟分析，该弛豫方案可作为产生相模拟的基础。

图14 水合物分解模拟中6.5 ns (a)及7.3 ns (b)时体系构型Fig. 14 System configuration at 6.5 ns (a) and 7.3 ns (b)during the hydrate decomposition simulation

3 结论

在进行产生相分子模拟前，需要对体系进行弛豫以合理化初始结构。对于不同体系与模拟需求，弛豫过程的设置不尽相同。

（1）需选取适宜的弛豫时间。一定时间弛豫后，构型变化并不明显，过长弛豫时间会显著增加计算时长。

（2）在NVT系综下，弛豫有利于提高体系的结晶度，更低的弛豫温度可更有效优化并稳定构型，但温度的设置需综合考虑所选力场模型的相平衡条件及冰点。在体系基本稳定后（合理的弛豫时间下），温度影响不再显著。故应结合后续模拟需求考虑，选择适中的弛豫温度。

（3）NPT系综下的弛豫有利于水合物结构的形成，越高的弛豫压力可在越短弛豫时间内得到可识别的笼型结构，但其稳定性越差。弛豫压力应根据后续模拟压力需求，结合相平衡条件及选定的弛豫温度而确定。

（4）不同弛豫系综所产生的弛豫效果不同。在稳定温度的基础上，NVT弛豫更利于结构的稳定，而NPT系综对笼型结构的形成具有促进作用。在体系弛豫时可优先考虑NVT弛豫以平衡温度与稳定结构，再结合NPT弛豫以促进水合物结构达到更接近理论数量的稳定状态。

因此，对于本文3 × 3 × 3 sI型CH4水合物晶胞，250 K的弛豫温度及10 MPa的弛豫压力可以得到较好的弛豫效果。在此温度和压力条件下，首先进行100 ps NVT弛豫，再进行100 ps NPT弛豫，即共200 ps的弛豫处理，所得到的平衡结构可基本拟合水合物分解温度及分解过程，具备开展产生相分子模拟的适用性。

符号说明

p―弛豫压力，MPa；T―弛豫温度，K；q―原子电荷值，e；σ―L-J势函数参数，nm；ε―L-J势函数参数，kJ/mol。