0 引言

分子动力学模拟是依靠牛顿力学来模拟分子体系运动的一套分子模拟方法，不仅可以用于模拟物质的宏观性质，还可以得出微观结构、粒子运动等信息^[1]。一般认为压强超过100 MPa就是超高压，物质在高压作用下，其物理和化学性质会发生巨大的变化，在超高压作用下，液体的压缩性显著增大^[2-3]。高压物理实验面临着实验条件、费用、实验设备的限制，而分子动力学计算利用计算机模拟，可以很好地解决这些问题。乙醇分子是常见的有机、极性分子，常用于化工生产中，研究其热力学性质和动力学性质对工业生产有重要意义。且乙醇和水在很多性质上都存在相似性，被认为是水的理想替代物^[4]，关于乙醇性质的研究也一直是热点。张阳^[5]用Monte Carlo方法研究了超临界乙醇的结构性质和氢键，压强范围从常压到高压70 MPa；Saiz等^[6]用分子动力学方法模拟研究了不同热力学状态下液态乙醇的结构、氢键和扩散等动态性质；Markus等^[7]研究了温度和压强(0.1~35 MPa)对甲醇和乙醇体系中氢键的影响。目前对乙醇的研究压强一般都在100 MPa以内，更高压强下乙醇性质的文章还少见报道。因此本文对液态乙醇的研究进行了补充，采用平衡分子动力学方法，模拟温度298~500 K、压强从10 MPa到超高压300 MPa条件下不同状态点乙醇的热力学性质、结构性质和动力学性质。本文的乙醇分子动力学模拟值可为其传质性质应用提供理论参考。

1 模型建立 1.1 势能模型

势能模型的选取对模拟结果的准确与否有着决定性影响。目前文献中分子模拟常用的乙醇模型主要是全原子模型和联合原子模型。鉴于全原子模型所需的计算时间较长，因此大部分文献采用联合原子模型进行模拟。联合原子模型有OPLS-UA(optimized potentials for liquid simulations-united atom)模型和TraPPE-UA(transferable potentials for phase equilibria-united atom)模型两种。OPLS-UA模型由Jorgensen^[8]开发，主要针对液相有机小分子，能够很好地模拟液态乙醇的各种动态性质。但是，在Chalaris和Samios^[9-10]用OPLS-UA模型对超临界甲醇进行分子动力学模拟的研究中，压强模拟值的准确性很差，说明OPLS-UA势能模型在压强较高时并不能真实反映分子的动态性质。TraPPE-UA模型是由Siepmann^[11]小组开发，由于其模型参数已经用实验的相平衡数据优化过，因此在很大的压强范围内都有较好的准确度。本文采用TraPPE-UA势能模型对高压下乙醇分子团簇的性质进行研究。

在TraPPE-UA模型中，乙醇分子由CH₃、CH₂基团和O原子、H原子合计4个作用位点构成。乙醇分子间的势能计算式为

式中：q_i、q_j——第i、j个原子或基团上所带电荷；

r_ij——原子或离子i与j间的距离；

σ_ij、ε_ij——原子间L-J作用参数。

式(1)右侧第1项表示短程Lennard-Jones(L-J)势，第2项表示长程库仑势。TraPPE-UA模型和OPLS-UA模型的相关参数见表 1。

1.2 实验方法

本研究采用Materials Explorer软件，在温度为298，400，500 K，压强分别为10，100，300 MPa条件下对乙醇分子体系进行了分子动力学模拟。采用Nose-Hoover热浴法控制温度，完成不同压强条件下的等温模拟^[13]。分子动力学模拟需进行以下假设：1)两体有效势能近似；2)最小镜像准则^[14]。分子动力学模拟的起始构型为面心立方晶格，将乙醇分子的起始取向定为随机分布，采用立方周期性边界条件。模拟过程中所用的乙醇分子总数为100个。采用球形截去法进行位能截断^[14]，采用Ewald加和法求解长程静电相互作用力^[15]，采用Shake算法固定乙醇分子的几何构型^[16]。分子动力学模拟过程中，每个分子都可视为是刚性体，运用五阶Gear预测-校正法求解体系的运动方程^[17]。分子动力学模拟过程中，时间步长取10^-15 s，单次模拟总时间取3×10^-10 s，采用NVT系综进行10^-10 s的初始化，使得体系达到相对平衡，随后在NPT系综中完成2×10^-10 s的模拟，完成各种参数的计算，模拟过程中每隔10^-14 s就输出一次构型，用于结果分析。

2 结果与分析 2.1 热力学性质

焓是物质重要的热力学状态量，体系焓计算式为

式中：E——体系的内能，kJ；

P——体系的压强，MPa；

V——体系的体积，m³。

表 2是本文模拟温度为298，400，500 K，压强为10，100，300 MPa时乙醇的焓值，并与文献[18]的模拟值和文献[19]用状态方程推导的焓值进行比较，结果较为接近，可见本文模拟结果可靠。

从表中可以看出，当压强相同时，随着温度的升高乙醇体系的焓值逐渐增大。这是由于当温度升高，乙醇分子间的热运动加剧，体系内能增大，焓也跟着增大。相同温度下，本文模拟的焓值随压强的增大而增大，但文献[18-19]中焓随压强的变化关系不明确，由于现有文献中尚未见到高压下乙醇焓值实验数据的相关报道，此计算结果有待进一步验证。

2.2 结构性质

径向分布函数(RDF)是反映物质微观结构的重要物理量，主要用于描述距离中心粒子为r处出现另一粒子的概率密度与随机分布概率密度的比值^[13]，其计算式为

式中：ρ——系统密度；

dN——距离中心分子(r，r+dr)区间内的分子数目。

图 1为T=298 K时乙醇分子的氧氢径向分布函数随压强的变化。图 2为P=10 MPa时氧氢径向分布函数随温度的变化。

从图 1可看出，压强变化对OH径向分布函数第1峰影响较大，随着压强的增大，OH径向分布函数第1峰分别出现在1.69，1.67，1.65Ǻ，峰位逐渐左移，峰值依次减小，第2峰也下降，但不如第1峰明显，在Petravic^[20]和Dellis^[21]等的研究中，也发现了OH径向分布函数的第1峰随压强增大而减小的现象。计算结果表明，随着压强的增加，第一配位圈内乙醇分子间的氢原子与氧原子间距逐渐减小，氢键作用逐渐增强。

从图 2可看出，当压强相同时，随着温度的升高，OH径向分布函数第1峰和第2峰逐渐下降，峰谷抬高，峰位右移。计算结果表明随着温度的升高，乙醇分子间的氢原子与氧原子间距逐渐增大，氢键作用逐渐减弱。表 3为各状态点乙醇的OH、OO、HH径向分布函数第1特征峰值与峰位。随着温度的升高和压强的增大，OH、OO、HH径向分布函数第1峰峰值均下降，但OH径向分布函数的变化较OO、HH径向分布函数更为明显。

图 3、图 4为100 MPa时OO、HH径向分布函数随温度的变化。从图中可看出，相同压强下，随着温度的升高，OO、HH径向分布函数第1峰下降，峰位右移，峰谷抬高，第2峰也下降，且在500 K时，第2峰消失。这说明随着温度的升高，乙醇的长程有序程度逐渐下降。这与高压下水的结构随温度和压强变化的趋势相同^[22]。

图 5~图 8为乙醇的MeH、MeO、MeMe径向分布函数。随着温度的升高和压强的增大，各径向分布函数第1峰变化最明显，峰高下降，峰谷抬高。500 K时，MeH径向分布函数的第2、3、4峰峰高已超过第1峰。MeMe径向分布函数第2峰逐渐消失。

通过对比不同原子和基团间径向分布函数的第一波谷位置r_min1可发现，r_OH＜r_HH＜r_OO＜r_MeH＜r_MeO＜r_MeMe。计算结果表明：相同温度和压强下，第一配位圈内亲水的O原子与H原子距离最小，而疏水的烷基之间的距离最远。

2.3 动力学性质

结合国内外研究现状，自扩散系数D的计算方法主要包括Einstein法和Green-Kubo法。Eintein法通过对均方位移(MSD)求斜率得到扩散系数。鉴于分子动力系统中的原子在永不停息地做无规则运动，各原子在不同时刻的位置都不相同。将(t)、(0)作为粒子在t时刻和初始时刻的位置，则粒子均方位移计算式^[13]为

式中 < >表示平均值。

根据统计学原理，只要分子数目足够多，计算时间足够长，系统的任一瞬间都可以当作时间的零点，所计算的平均值应该相同。根据爱因斯坦的扩散定律，扩散系数(diffusion constant)与均方位移的关系是：

式中D即为粒子的扩散系数。

当模拟的时间很长时，均方位移对时间曲线的斜率的1/6就是扩散系数。表 4是本文用均方位移计算的乙醇的自扩散系数。

由表可知，在相同压强下，随温度的升高，乙醇的自扩散系数增大。这与文献[23]测定的常压下乙醇自扩散系数随温度的变化趋势相同，且随着压强的增大，这种变化趋势逐渐减缓。而相同温度下，随着压强的增加，乙醇自扩散系数逐渐减小。这是因为压强增大使乙醇分子之间的距离变小，分子间排列更为紧密，相互碰撞现象明显增加，从而阻碍了分子的扩散运动；当温度不断上升时，乙醇分子间的距离逐渐增加，相互作用力减弱，此时将有利于分子的扩散运动。本文模拟值与文献[18]的模拟结果基本一致，但在500 K，10 MPa时相差较大，由于实验条件的限制，高压下乙醇分子团簇的自扩散系数实验数据匮乏，有待进一步验证。

3 结束语

本文模拟了温度298~500 K，压强从高压10 MPa到超高压300 MPa条件下乙醇分子体系的热力学性质、结构性质和动力学性质，并得出以下结论：

1) 相同压强下，乙醇体系的焓值随温度的升高而增大，本文模拟值与文献模拟值接近，本文所用TraPPE-UA模型可靠。

2) 通过对径向分布函数的研究，发现在压强很大的情况下，乙醇体系的结构依然很有规律。随着压强和温度的升高，各径向分布函数的第1峰高度明显下降，乙醇的长程有序程度下降。

3) 当压强相同时，乙醇的自扩散系数随温度的升高而增大，但压强增大时，这种变化趋势逐渐减缓。乙醇的自扩散系数随温度的升高而增大，而压强对扩散系数的影响与温度正好相反。