超临界水中氢键特性及其空间分布的密度相关性

2020-08-07曾冬雷宋嘉文范利武

高校化学工程学报 2020年3期

曾冬雷, 冯飙, 宋嘉文, 范利武,2

(1. 浙江大学热工与动力系统研究所, 浙江杭州 310027;2. 浙江大学能源清洁利用国家重点实验室, 浙江杭州 310027)

1 前言

超临界流体是温度高于临界温度Tc与压力高于其临界压力pc的特殊流体，由于其普遍具有较高的氧化性、溶解度和扩散系数，因此常被用做氧化剂[1-2]或催化剂[3]等反应介质。近年来，随着对超临界流体物性研究的不断深入，其密度的局部不均匀性受到了研究人员的广泛关注，并采用分子动力学方法进行了研究。分子动力学方法可从微观机理角度解释宏观现象并且不受实验技术限制，可模拟极端条件下的体系并计算得到相关物理量[4]，是研究以超临界水(临界温度与压力分别为647.15 K 与22.1 MPa)为代表的超临界流体物性的有力工具。

SKARMOUTSOS 等[5]发现温度为 666 K 的超临界水(模拟的水分子数量为 500 个，密度范围为 0.2 ρc～2 ρc，ρc为水的临界密度)存在局部密度与结构的不均匀性，而且含氢键流体的局部密度不均匀性更大[5-6]。氢键的键能介于共价键和范德华力之间，考虑分子间氢键的平均强度，形成一个氢键则增强了其周围氢键的协同作用，且氢键的协同作用可使水分子以大的团簇结构存在[7]。MA 等[8]从氢键数量与体系能量波动来解释这种不均匀性，并发现温度在666 K 的条件下，密度在低于0.19 ρc时能量波动急剧增加，此时四面体序参数有最小值，体系短程有序被破坏。SKARMOUTSOS 等[9]发现当温度为666 K 时，周围有1 个和2 个氢键的水分子对超临界水的局部密度不均匀性的贡献最大，并从三角形和四面体序参数角度揭示了局部有序性与氢键的关系。SAHLE 等[10]则通过从头分子动力学模拟和密度泛函理论方法，发现在温度低于873 K 与压力低于132 MPa 的超临界水中仍然存在着氢键网络。

由于密度不均匀性可能引起超临界水中氢键特性的变化，并且氢键作为贡献较大的分子间相互作用，其数量、键能大小及空间分布会对超临界水的物性造成重要影响。因此对氢键特性及其空间分布随密度变化的研究尤为重要。但是目前文献中关于超临界水密度对氢键键能及其空间分布影响的研究较为缺乏。因此，本文参照现有文献的研究结果，对温度为近临界温度666 K 下不同密度的超临界水体系进行分子动力学模拟，并增加更高温度的700 与800 K 这2 个工况作为对照组，研究不同温度下超临界水的密度对其氢键数量、键能以及空间分布规律的影响，将有助于进一步理解超临界水密度不均匀性的微观机制。

2 分子动力学模拟

2.1 模拟细节

本研究采用经典的TIP4P/2005[11]作为水分子模型。如图1 所示为超临界水的初始分子建模示意图，其中红球表示氧原子，绿球表示氢原子，立方体的边长为L。所模拟的超临界水体系含1 000 个水分子，不同密度下模型的L 取值如表1 所示(ρ/ρc为超临界水密度与水临界密度的比值)。

分子动力学模拟使用的平台是LAMMPS，采用的力场为TIP4P 力场，并设置三维周期性边界条件(PPP)。温度与压力的控制采用Nose-Hoover 耦合方法，并采用PPPM 算法计算静电相互作用力。初始分子速度由高斯随机分布决定。库伦截断半径与范德华力截断半径均为1 nm。根据相关文献中对水分子进行模拟的经验，模拟的时间步长设为10-15s，每隔10-12s 保存坐标数据。在正则系综(即NVT 系综)下进行模拟，系统先运行2×10-9s 达到平衡态，再计算10-9s 进行数据分析。如上所述，为揭示不同温度的影响，系统的温度设定为666、700、800 K。

图1 超临界水的初始分子建模Fig.1 Initial molecular model of supercritical water

表1 超临界水模型的边长取值Table 1 Side length values of the supercritical water model

2.2 模型检验

为了对模拟体系的分子数量无关性进行检验，首先计算了温度为666 K 时分子数分别为500、1 000与5 000 的3 个超临界水体系中的平均氢键数随密度的变化关系。如图2 所示，3 个体系模拟得到的平均氢键数(＜ nHB＞permolecule)随密度比值的变化规律完全一致，体系分子数的变化对模拟结果影响不大，相对偏差在20%以内。因此，为了相对节约计算资源，本研究的后续模拟中均选取水分子数为1 000 来构建模拟体系。

对模拟得到的超临界水的自扩散系数与文献中的实验值[12]以及前人的分子动力学模拟结果(采用 256个水分子，温度为666 K)[8]进行对照，以检验本文所建立的超临界水模型的可靠性。采用Einstein 法[13]计算超临界水的自扩散系数。该方法的基本原理是通过统计粒子均方根位移与运动时间的关系，当时间足够长二者呈线性关系时取其斜率的1/6 即为自扩散系数，其公式如下所示：

图2 不同水分子数超临界水体系中平均氢键数随密度的变化Fig.2 Variation of average number of hydrogen bonds as a function of density in supercritical water systems with different water molecule numbers

式中：N 为粒子数目，t 为运动时间，ri(t)为粒子瞬时位置，ri(0)为粒子初始位置。

实验值由LAMB 等[12]的实验数据插值计算得到。如图3 所示，本文的模拟值与实验值和前人的模拟值均吻合较好，总体偏差率低于15%。在保持温度不变的条件下，超临界水的自扩散系数D 随密度比值升高而减小至趋于平缓。这是由于密度增大，分子间的碰撞概率增大，且氢键数量增加，进一步阻碍了分子扩散，传质能力减弱，从而表现为自扩散系数减小。在ρ ＞ 1.2 ρc时自扩散系数的变化低于 10%，与 MA 等[8]在ρ ＞ 1.25 ρc时发现的密度对自扩散系数影响很小的结果基本一致。800 K 超临界水自扩散系数最高，666 K 最低，且随着温度的降低自扩散系数也逐渐减小，是由于温度升高促进了水分子的热运动，且氢键作用减弱，水分子运动的空间位阻效应减弱[14]，传质能力增强。

2.3 径向分布函数

如图4 所示为径向分布示意图，通常用来描述局部粒子的分布，其表达式如下所示：

图3 超临界水自扩散系数随密度的变化Fig.3 Variation of self-diffusion coefficient of supercritical water as a function of density

图4 径向分布示意图Fig.4 Schematic diagram of radial distribution

式中：dnr为半径从r 到r+dr 球壳内围绕每个A 粒子周围B 粒子的数目，ρB为B 粒子的数密度，r 为距中心原子的半径，dr 为球壳厚度。径向分布函数实际上是A 粒子周围B 粒子出现的概率密度函数，可通过其观测粒子的空间分布。由于每个水分子中只包含 1 个氧原子，因此可以通过观察氧氧径向分布函数来描述水分子的分布状况，即从 gOO(r)峰值出现的位置，高度及曲线的走势来分析水分子的分布状况。因为水分子的分布状况的变化会影响氢键的分布状况，所以分析水分子分布对氢键分布的研究具有重要意义。

2.4 氢键特性表征

本文使用几何标准定义氢键[6]，从而对计算输出的坐标文件进行编程统计氢键的数量。如图5 所示，供体与受体氧原子距离rOO≤0.36 nm，供体与受体氧原子与形成氢键的氢原子夹角∠HOO≤30°且受体氧原子与氢原子距离rOH≤0.24 nm 可视为形成氢键。

图5 氢键几何模型Fig.5 Geometric model of hydrogen bonds

氢键键能仅由3 个原子(OD-H...OA)的短程库仑作用决定：

3 模拟结果及讨论

3.1 水分子的空间分布

由于每个水分子只含有一个氧原子，因此可以通过计算氧氧径向分布来分析水分子的空间分布情况。不同密度超临界水中的氧氧径向分布模拟结果如图6 所示。

图6 超临界水中的氧氧径向分布Fig.6 Radial distribution of oxygen-oxygen in supercritical water

根据NVT 系综下分子运动轨迹文件，计算了在666，700与800 K 温度下的不同密度的氧氧原子的径向分布函数。由图6可知，超临界水gOO(r)的第1 个峰均位于0.29 nm，说明超临界条件下水分子之间出现概率最高为0.29 nm。第2 峰基本消失，说明在超临界条件下氢键作用变得很弱。在恒温条件下，gOO(r)的第1 峰高度具有密度依赖性，且随着密度的减小，第1 峰高度增大，说明超临界水短程有序性随密度减小而增大，在较低密度的条件下会形成局部团聚现象。水分子的局部团聚会导致氢键在空间分布上呈现局部团聚的现象。但随着温度升高，由于分子热运动增强，这种密度依赖性减弱，且第1 峰高度减小，这与SOPER 等[15]的研究结论一致。

3.2 氢键数量与平均键能

水分子平均氢键数(＜ nHB＞permolecule)与密度比值的关系如图7 所示。值得注意的是，图中给出了在若干次模拟中所得到数据点的误差波动情况以说明模拟数据的可重复性。由于使用了NVT 系综，因此在模拟过程中存在体系温度的随机波动现象，从而导致模拟结果也会出现一定的随机波动。在温度不变的情况下，水分子平均氢键数随密度增加而增加，这是由于水分子密度增大，结构更紧凑，水分子之间更容易形成氢键。但当 ρ ＞1.4ρc时平均氢键数趋于平缓，密度的影响减小。超临界水在666 K 下平均氢键数最少，800 K时最多。这是由于水分子的热运动增强导致更难达到成键的几何条件。因此，随着温度的升高，氢键的作用减弱。在666 K 温度下，平均氢键数在密度为2.0ρc时的平均氢键数较高，可能是因为在较低温度下氢键的密度依赖性更大，高密度时更容易出现团聚现象。

图7 超临界水中氢键数量随密度的变化Fig.7 Variation of hydrogen bond number as a function of density in supercritical water

图8 超临界水中平均氢键键能随密度的变化Fig.8 Variation of average hydrogen bond energy as a function of density in supercritical water

平均氢键键能EHB与密度比值的关系如图8 所示。在温度一定的条件下，超临界水平均氢键键能随着密度比值增大有上升的趋势，其中温度为666 K 时平均氢键键能在2.0 ρc相比0.4 ρc增加了7.0%，温度为700 K 时增加了9.1%，温度为800 K 时增加了12.3%。氢键键能增长幅度随温度升高而增加，总体增长幅度较为平缓。这是因为密度增大导致水分子之间距离更近，形成的氢键更稳定，键能更大。在温度为800 K 的条件下，超临界水平均氢键键能则波动最大，可能是由于水分子热运动较强，使得原子之间的距离有很大的随机性；在温度为666 K 的条件下，平均氢键键能总体最高，且随和温度的上升而下降，这是因为温度更低的条件下分子间距离更小，导致氢键键能更大。

3.3 氢键的空间分布

图9 不同密度(0.4 ρc，1.2 ρc，2.0 ρc)超临界水中氢键的空间分布Fig.9 Spatial distribution of hydrogen bonds in supercritical water having different densities (0.4 ρc，1.2 ρc，2.0 ρc)

如图9 所示以受体氧原子坐标为氢键坐标，氢键键能大小表现在圆球的大小上。温度不变时，随着密度的增大，氢键成键在空间分布上更均匀，团聚现象(图中黑色虚线所示)越少，密度增大与氢键更均匀的分布使得水分子的传质能力降低，表现在自扩散系数的减小。高密度超临界水体系中心处平均氢键键能大小更高(圆球更大)，成键更稳定。在密度不变的条件下，氢键空间分布随着温度的升高而变得更加均匀，即更难出现团聚现象。同时，在更高温度的情况下，氢键数目更低，氢键键能大小差异更小。