超临界Lennard-Jones流体结构特性分子动力学研究*

王艳 徐进良2)† 李文 刘欢

1) (华北电力大学, 低品位能源多相流与传热北京市重点实验室, 北京 102206)

2) (华北电力大学, 电站能量传递转化与系统教育部重点实验室, 北京 102206)

研究超临界流体在不同压力和温度的结构特征有助于深刻理解并有效利用超临界流体.本文采用分子动力学方法模拟超临界压力、拟临界温度附近流体的结构及密度波动曲线的排列熵, 分析状态参数变化的影响.结果表明, 定压下, 径向分布函数随温度升高, 第一峰值位置逐渐向右移动, 但右移幅度随着压力偏离临界点距离的增大而减弱, 近临界压力时, 出现峰值最高点的工况和等温压缩系数的极值点位置一致, 压力增大, 该现象消失.低压力拟临界点时易出现面积大、相对集中且分布稳定的高/低密度区, 无明显嵌套现象.静态结构因子存在一定发散行为, 发散的最大值和等温压缩系数极值点所处工况符合.低压力时密度时间序列的波动幅度最大, 类周期现象较明显.在分子间势能、等温压缩系数和热运动效应的共同作用下, 当压力(P)为 1.1 倍的临界压力 (Pc)时, 排列熵在 0.99 倍的拟临界温度 (Tpc)达到最小值, P = 1.3Pc和 1.5Pc时, 最小排列熵与等温压缩系数的最大值工况点保持一致, 压力继续增大, 各模拟工况密度和排列熵的波动减弱,流体均匀性增强.

1 引 言

严格意义上超临界流体是指温度和压力均在临界点以上的流体[1,2], 传统教科书认为超临界态是一个均匀状态, 没有气液区分.然而, 随着超临界流体在工业中广泛的应用, 越来越多的研究得到开展.已有研究[3,4]发现, 可以根据定压比热容极值点的位置将超临界区分为类液和类气区, 即Widom线, 利用速度自相关函数的突变点的连线,得到的Frenkel线是类液类气区分界线的另一个划分标准[5,6].此外, Banuti等[7,8]采用分子动力学和理论推导的方法进一步将超临界流体划分为类液、类气和过渡三个区域, 而分子动力学模拟的方法以其成本低、安全、便捷等优点在超临界流体模拟中得到广泛的应用.

目前关于超临界流体的分子动力学模拟主要集中在热物性和输运性质的计算以及不同势函数、截断半径选取对计算结果的影响等方面[9−12], 为后续相关研究提供可靠的设置参数.Skarmoutsos和Samios[13,14]采用分子动力学方法研究超临界水沿等温线, 系统温度( T )和临界温度( Tc)的比值T/Tc=1.03, 局部密度增强和动力学特性与系统密度的依赖关系, 此外还对该等温线上甲醇和CO2的密度增强效应进行研究, 发现两种流体均在0.7倍的临界密度( ρc)时得到最大的增强效应;Yoshii和 Okazaki[15]对 Lennard-Jones (LJ) 流体沿等温线 T /Tc=1.03 , 多个密度下超临界流体的结构和团簇进行研究, 认为静态结构因子在临界区域内存在较强的发散行为, 临界密度附近存在临界慢化现象[16]; Metatla等[17]根据径向分布函数、计算配位数与期望配位数之间的关系研究确定400 ℃超临界水是存在高/低密度区的非均质结构;Maddox等[18]采用分子动力学方法模拟了二维LJ流体在接近临界点时的密度不均匀特性, 指出该现象主要是“势能诱导”和“临界波动”两种作用机制共同影响的结果; Yamane等[19]不但分析了径向分布函数和静态结构因子分布曲线与温度的关系,而且研究了近临界点附近流体汞的配位数和截断半径的依赖关系, 得到在近临界点附近流体结构分子动力学模拟需要使用大规模的系统和较大的截断半径.除分子动力学外还有部分学者采用拉曼散射和小角度X射线散射的实验方法对超临界流体的非均匀性进行大量的研究[20−23].

综上所述, 现有文献对超临界流体的结构特性系统的研究较少, 文献几乎都集中在单一的等温线上, 涉及到状态点较少, 在相图上的位置不明确,而且没有考虑在跨越Widom线前后的结构特性的转变.即在超临界压力下, 拟临界点附近流体的结构特性和随时间演化过程的研究仍处于空白.因此, 针对超临界压力拟临界点附近超临界流体的结构和时变特性的研究是非常必要的.本文利用分子动力学方法模拟研究不同压力、不同温度下超临界流体的结构特点及空间和时间演化特性, 分析压力、温度对相关物理参量的影响, 首次在分子动力学模拟中引用排列熵的概念对密度时间序列曲线进行分析, 剖析不同工况下高/低密度区形成及产生最小排列熵的主要作用机制.研究结果为从微尺度层面揭示超临界流体的特性提供可靠支撑, 也对超临界流体的实际应用提供有益启发.本模拟采用开源分子动力学模拟软件LAMMPS 实现, 分子位型采用Ovito软件进行可视化.

2 物理模型及模拟细节

图1(a)为物理模型示意图.系统沿着 x, y,z三个方向均采用周期性边界条件.模拟体系尺寸为 Lx=Ly=Lz=L , 模拟系统内充满超临界流体氩, 初始氩原子按照FCC晶格排列方式布置.有文献[19]研究表明, 超临界流体模拟过程中需要较大的系统尺寸, 得到的计算结果具有较高的可靠性, 因此本文模拟过程中, 各模拟工况均维持系统原子数为104个, 根据不同计算工况的温度和压力,确定各工况对应的密度, 得到模拟体系的尺寸范围为27.2283 σ —40.8040 σ.

超临界流体氩的临界点温度(Tc)为150.687 K,压力 ( Pc)为4.863 MPa, 密度 ( ρc)为 535.6 kg/m3,为揭示不同的超压力、拟临界温度(Tpc)附近超临界流体结构及密度时间序列曲线波动特性, 选择1.1 Pc—2.0 Pc4 个超临界压力, 如图1(b)所示, 每个压力下取0.95 Tpc—1.05 Tpc7个工况进行模拟分析.

图1 (a) 物理模型图; (b) 模拟状态点在相图中的分布Fig.1.(a) Physical model of system; (b) simulation points on phase diagram with Widom line, liquid-like and gas-like region.

图2 (a) 定压比热容 (cp)变化曲线; (b) 等温压缩系数 (kT)变化曲线Fig.2.(a) The curve of cp under different pressures; (b) the curve of kT under different pressures.

根据超临界流体的物性参数和相关研究结果可以得到, 定压比热容(cp)的极值点连成的线即为Widom线, 也是超临界区类液类气的分界曲线,计算工况在定压比热容的曲线上的分布如图2(a)所示, 计算工况包括了类液和类气区域, 且都集中在拟临界点附近.相应的等温压缩系数(kT)的变化曲线如图2(b)所示, 从图中可以得到在P=1.1Pc时, 等温压缩系数在拟临界温度下取得极值, 随着压力的增大, 取得极值点的位置发生右偏, 在P=1.3Pc和 1.5Pc时, 均在 1.01Tpc得到极值点, 当压力增大到 2.0Pc时, 等温压缩系数在文中选取的计算温度区间内呈现一个平缓的上升趋势, 没有出现极值点.超临界流体物性突变是流体结构发生变化的具体表现, 也是各工况计算分析中需要重点考虑的因素.

分子动力学的基础是牛顿第二定律, 直接给出了加速度和施加外力之间的关系:

其中m和 r 分别为分子的质量和矢量位置, 下标i表示原子i, 右端的二阶位移导数项表示原子i的加速度 ai, Fi为原子i所受的总作用力.

流体氩分子之间的相互作用采用Lennard-Jones (LJ) 势能模型, 表达式为

其中r为原子对间的距离, 液体氩原子之间的尺寸参数 σ =0.3405 nm , 能量参数 ε = 1.67 × 10–21J,分子质量 m = 6.69 × 10–23g.根据相关文献研究,在超临界模拟过程中势能截断半径为5.88 σ[24], 超过此距离的分子, 其相互作用忽略不计.

采用Velocity-Verlet算法求解运动方程, 时间步长取为1 fs.在模拟过程中对整个系统施加NVT 系综, 采用 Nose-Hoover控温方法.每个算例计算 15 ns, 前 10 ns为充分弛豫平衡过程, 后5 ns用于统计分析各种参量, 后续的统计过程中,忽略弛豫过程的时间, 直接从统计时刻开始记录时间为 0τ.

3 结果与讨论

3.1 径向分布函数及配位数

径向分布函数(RDF)是系统的区域密度与平均密度之比, 分子动力学中计算径向分布函数的方法为[25]

其中 ρ 为系统的密度, N为分子的总数目, D为计算的总时间 (步数), δr 为设定的距离差, 参考分子与其中心的距离 rc→ rc+δr 间的分子数目为 ∆ N.

图3 径向分布函数 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0PcFig.3.Radial distribution function: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.

由图3得知, 在各计算工况下, 温度低于拟临界点温度时, 径向分布函数会出现高低不同的峰谷, 在1附近产生振荡, 表现为一种“短程有序、长程无序”的类液体现象, 和文献[26]中描述的液相趋势一致.随着温度的升高, 波谷逐渐抬升, 第二波峰逐渐变小, 在 T =1.05Tpc时第二个峰值基本消失, 满足气相在短程区域有一个峰值, 而后单调下降到1的变化趋势, 在各模拟工况下, 随温度的升高径向分布函数表现为类液-类气的特征, 两种区域的结构存在差异.在定压下, 随着温度的升高,密度逐渐降低, 第一峰值出现的位置逐渐向右偏移, 随压力的增大, 偏移程度逐渐减小, 在足够大的压力下该偏移量将会消失.四个分图对比发现,随着压力的增大, 第一峰值逐渐减小, 峰宽度也减小, 这主要是由于超临界流体氩局部出现聚集的高密度区的结果.氩的局部聚集必然会导致体系内部出现“空隙”, 压力的增大将这些“空隙”压缩, 减小了局部聚集氩相对密度, 因而会出现峰值和峰宽度都减小的现象, 和现有文献[27]分析保持一致.同时也进一步证明超临界流体的高压缩性.在P=1.1Pc和 1.3Pc时, 随着温度的升高峰值呈现出先下降、而后升高再下降的一种变化规律, 在T=1.0Tpc和 1.01Tpc时得到最大值; 在 P =1.5Pc和2.0Pc时, 随着温度的增加, 峰值整体呈现一个下降趋势, 在 T =1.01Tpc时仍呈现微弱的增大趋势, 但是峰值最大值仍出现在最低温度工况下, 这主要由于在拟临界点附近超临界流体存在很大的压缩性,促使流体聚集产生高密度区, 但是随着压力的增加, 流体的压缩性迅速降低, 温度升高带来的分子热运动逐渐加强, 使得流体的聚集能力逐渐减弱.

配位数(CN)的变化趋势和密度呈线性关系,也是描述流体微观结构的重要物理参数, 反映了距离中心分子为rx的球体区域内分子的个数, 描述了分子排列的紧密程度, 配位数越大, 粒子排列越紧密, 不同区间范围内的配位数的定义计算式为[28]

计算超临界流体氩的配位数能清楚反映出流体的内部结构.具体模拟结果如图4所示.

有研究表明配位数不仅与结构、密度有关, 而且是温度的函数[29].在相同的压力下, 随着温度的升高, 流体的密度减小, 此时分子间的距离增大,分子间引力对配位数的影响越来越大, 但是温度的升高, 促使分子的热运动加剧, 在引力和热运动的共同作用下, 发现随温度增加配位数逐渐减小.此外随着压力的增大, 确定的温度区间内, 配位数的波动范围逐渐缩小, 各工况的差异性逐渐减弱.从图5 可以观察到, P =1.1Pc时, 随温度的增加, 配位数曲线的斜率在16.22—5.90区间变化, 随着压力的升高; P =2.0Pc时, 曲线斜率在 14.27—9.25区间变化.压力越接近临界点时, 随着温度的变化,密度波动范围越大, 从而引起配位数的波动区间较大, 而在远离临界点的压力下, 密度随温度的变化范围较小, 相应的配位数曲线的斜率范围较窄.在T/Tpc<1.0的类液区, 压力增大, 拟临界点温度升高, 密度减小, 配位数斜率呈现减小的变化趋势;在 T /Tpc≥1.0 的类气区表现为相反的变化趋势.可以得到密度是影响流体配位数的主要参数, 温度和压力也是通过调整密度的大小间接影响配位数分布.

图4 配位数 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0PcFig.4.Coordination number: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.

根据配位数的计算值( Nc)和期望值( Ne)对密度不均匀性进行判断, 如果直接比较两个数字的大小, 这样的标准过于严格, 局部结构的微小波动,即使只有一个分子穿过设置的半径边界, 也会导致密度的分类从平均转变成高或低密度.根据文献[30]提到的稍微宽松且能准确判断的标准, 即选取一个小量的值 δ , 用 Ne± δ 作为参数度量平均密度的变化, 具体划分原则如下:

图5 配位数曲线斜率Fig.5.The slope of coordination number curve.

计算中允许的波动量为30% Ne, 则 δ 应满足2 δ +1=0.3 Ne, 进一步得到不同参数下 δ 的具体值,利用划分原则可以判断密度分布趋势, 称该方法为“30% 方法”.采用该方法, 根据配位数得到在P=1.1Pc和 2.0Pc, T =Tpc流体在 xy 平面内高/低及平均密度区分布随时间的演化如图6所示.

图6(a)表示在 0—5000τ的时间范围内P=1.1Pc, T =Tpc时的演化过程, 由图可知, 该计算工况产生的均值区面积较小, 低密度区的位置主要集中在模拟系统的中部, 随着时间的演化, 密度在不停地波动, 低密度区呈现分散-聚合-分散的演化规律, 产生的高密度区的面积相对较大且位置集中, 演化过程中波动微弱.由于在超临界区表面张力的消失, 不存在亚临界工况下的弯曲界面.形成该现象的原因主要是由于当密度低时, 分子间引力起主导作用, 在温度的影响下分子处于不断的相互碰撞中, 能量的增加, 导致高密度区的形成.此外,该工况拟临界点位置出现等温压缩系数的极值点,形成较高的密度区所付出的代价变小, 超临界流体的关联长度在拟临界处也存在极大值, 因此会出现大面积的, 相对稳定的高密度区.在确定尺寸和分子数的系统中, 高密度区的形成必然会引起低密度区的产生.随着高密度区的形成, 分子间的距离缩短, 增大的斥力将部分分子从高密度区向低密度区推, 当分子间距离较大时, 引力起主导作用, 会使得分子再次聚集成高密度区, 但此时作用势效应相对等温压缩系数效应较小, 因此各区域所处位置稳定, 波动微弱.随着压力的升高, 流体的等温压缩系数仍存在极值点, 但是数值较小, 流体可压缩性减弱, 较难形成高密度区, 仅在分子间短程势作用下, 形成由几个分子组成的且较为分散的高密度区, 形成相对较大的平均值区, 如图6(b)所示.随时间的演化, 各区域不停波动, 存在明显的嵌套现象, 系统中局部出现类似“花斑”的现象, 这些花斑若隐若现、此起彼伏、互相嵌套的性质和相关教课书[31]中提出的结论保持一致.

图6 流体在 xy 平面内高/低密度区分布随时间的演化 (a) P = 1.1Pc, T = Tpc; (b) P = 2.0Pc, T = TpcFig.6.Liqud atoms evolution over the xy plane with different pressure: (a) P = 1.1Pc, T = Tpc; (b) P = 2.0Pc, T = Tpc.

图7 不同压力下, 拟临界点温度下高/低密度区占比 (a) 高密度区占比; (b) 低密度区占比Fig.7.The ratio of high/low density region at pseudo-critical point temperature under different pressure: (a) The ratio of high density region; (b) the ratio of low density region.

从图7可以直观地观察各工况不同密度区所占比例随时间的变化.不同的压力条件下, 类液和类气区的占比一直处于一个波动的状态, 仅观察曲线, 发现波动过程并没有实际的规律可循, 整体表现为一个混乱的动态变化过程.从图7(a)和图7(b)分图中发现随着压力增大类液和类气区的占比整体呈现下降的趋势, 进一步说明随压力增大, 形成高密度区较为困难, 计算压力距临界压力越远, 流体的均匀性越强.在 P =1.1Pc, T =Tpc的工况时,高密度区占比约为60%, 低密度区的占比约为35%,此时均匀区域的占比最小, 系统的不均匀性最强.

3.2 静态结构因子

结构因子主要表征材料对射线的散射能量, 反映材料结构的平均信息, 可以进一步应用到流体中, 观察流体的结构特性.而静态结构因子和径向分布函数互为傅里叶变换[19,32], 计算式为

其中 n = N/V, V 为系统体积, 由于结构因子是倒易空间, 变量k与距中心分子的距离 rc成反比.

由图8可知, 各模拟工况下均在较小的k值范围内存在发散行为, k < 0.5 σ−1, 随着压力的增加,这种发散行为逐渐减弱.在 P =1.1Pc时, 拟临界点处的发散行为最为强烈, 流体表现为较强的小角度散射, 曲线在一定范围内存在微小的波动, 而随着压力的升高, 这种波动现象消失.在 P =1.3Pc和1.5Pc时, 发散最强烈的行为出现在偏离拟临界点的 T =1.01Tpc工况, 与等温压缩系数的最大值点保持一致.在图8(a)—图8(c)分图中得到随着温度偏离拟临界点的距离增大, 发散性减小, 即各工况均在 T =0.95Tpc和 1.05Tpc时得到静态结构因子的最小值.在图8(d)分中, 由于等温压缩系数在一个较小的水平缓慢增加, 因此静态结构因子的发散特性也呈现微弱增加趋势.压力较低时, 静态结构因子的发散在拟临界点温度达到最大值, 但随着压力的增加, 发散的极值点工况与等温压缩系数的极值点工况相符合.

3.3 密度时间序列曲线和排列熵

密度随时间的变化是分子动力学模拟中比较容易得到的数据, 该曲线给出初步粗略的波动信息.从图9(a)可以看出, 压力 P =1.1Pc, 密度时间序列曲线呈现出包含类周期特征的大幅波动, 流体结构具有局部有序特征, 但这种波动存在不规则运动的相互叠加, 大波套小波的现象, 并不具有严格的周期性.随着压力的增加, 波动的幅度和类周期性均减弱, 随机性增强, 此时流体结构具有较强的无序性.

图8 静态结构因子 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0PcFig.8.Static structure fator: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.

图9 (a) 密度时间序列曲线 ( T = Tpc); (b) 排列熵Fig.9.(a) Time series of density for T = Tpc; (b) permutation entropy.

为具体描述密度时间序列的复杂程度, 引入排列熵的分析方法, 对一个确定长度的时间序列, 可以根据自相关函数法确定一个延迟时间 τD, 根据排列熵嵌入维度的要求, 文中选取各工况嵌入维数m = 5.排列熵作为衡量时间序列曲线复杂度的指标, 越规则的时间序列, 对应的排列熵越小, 越复杂的时间序列, 对应的排列熵越大.

根据文献[33]提出的计算方法, 对各工况的密度时间序列曲线进行排列熵的计算, 从图9(b)中可以观察到在 P =1.1Pc时, 在 T =0.99Tpc和1.0Tpc时排列熵的值较小, 对应较规则的时间序列.随着压力的增大, P =1.3Pc和 1.5Pc时, 均在T=1.01Tpc时得到最小的排列熵, 随着压力的进一步增加, 排列熵的值呈微弱增大趋势, 且各工况的值在一个水平线附近波动.产生这种现象的原因主要是分子间作用势、等温压缩系数和分子热运动效应共同作用的结果.P =1.1Pc时, 各计算工况的温度低, 分子热运动效应较弱, 等温压缩系数在拟临界点处达到极大值, 流体容易被压缩形成高密度区, 在这两种效应的共同作用下, 排列熵的最小值出现在T=0.99Tpc的工况.随着压力的升高( P =1.3Pc和1.5Pc), 拟临界温度逐渐增大, 分子的热运动加剧,导致分子间的作用势效应减小, 此时等温压缩系数效应在密度波动中占据主导地位, 最小排列熵和等温压缩系数极大值点工况一致.当 P =2.0Pc时, 各工况得到的排列熵变化趋势和等温压缩系数保持一致.

4 结 论

采用分子动力学方法模拟不同超临界压力, 拟临界温度附近流体的结构特征, 对流体的径向分布函数、配位数、静态结构因子、密度时间序列曲线及排列熵展开研究, 分析压力和温度的变化对各参量的影响, 得到如下结论:

1) 径向分布函数在类液区存在“短程有序、长程无序”的现象, 随着温度的升高, 这种有序的现象逐渐消失, 在类气区仅在短程区域存在一个峰值, 而后单调下降到; 定压下, 随着温度的升高, 第一峰值出现的位置逐渐右移, 但这种变化趋势随着压力偏离临界点距离的增大而减弱; 此外首次提出, 在 P =1.1Pc和 1.3Pc时, 峰值的最大值均在等温压缩极大值点工况出现, 但是随着压力的增加,该现象不再明显.

2) 定压时, 随着温度的升高, 配位数逐渐减小, 压力的增加导致不同温度下流体密度的差值减小, 因此引起配位数的波动范围缩小, 分布区域变窄; 采用配位数标定高/低密度区的标准, 在P=1.1Pc, T =Tpc时可以得到大面积、分布集中且波动较小的高/低密度区, 此时均值区域占比较小; 随着压力升高, 均值区域逐渐增加, 高/低密度区逐渐减小, 仅有几个分子大小, 且相互嵌套, 并随时间有较明显的波动, 说明随着压力的增大, 流体的均匀性逐渐增强.

3) 静态结构因子是通过散射效应观察流体的内部结构, 各模拟工况均在较小的k值内存在发散行为, 随着压力的增加, 静态结构因子曲线的发散值呈迅速下降的趋势, 研究发现定压时静态结构因子最大值和等温压缩系数极值点工况保持一致.

4) 密度时间序列曲线可以初步得到中心切片密度随压力的增加, 波动幅度逐渐减弱的规律, 和前文压力增大流体均匀性增强的结论符合.排列熵的变化规律可以分为三类: 一是在 P =1.1Pc的低压下, 排列熵在小于拟临界温度的位置得到最大值;二是在 P =1.3Pc和 1.5Pc的中压下, 排列熵的最小值点和等温压缩系数极值点一致; 三是在P =2.0Pc的高压下, 排列熵的值在一个水平线附近波动, 整体变化趋势较为平缓, 流体均匀性增强.探讨微尺度下超临界流体的结构特征有助于全面了解超临界流体特性, 为超临界流体的应用提供有力支撑.