-
研究超临界流体在不同压力和温度的结构特征有助于深刻理解并有效利用超临界流体. 本文采用分子动力学方法模拟超临界压力、拟临界温度附近流体的结构及密度波动曲线的排列熵, 分析状态参数变化的影响. 结果表明, 定压下, 径向分布函数随温度升高, 第一峰值位置逐渐向右移动, 但右移幅度随着压力偏离临界点距离的增大而减弱, 近临界压力时, 出现峰值最高点的工况和等温压缩系数的极值点位置一致, 压力增大, 该现象消失. 低压力拟临界点时易出现面积大、相对集中且分布稳定的高/低密度区, 无明显嵌套现象. 静态结构因子存在一定发散行为, 发散的最大值和等温压缩系数极值点所处工况符合. 低压力时密度时间序列的波动幅度最大, 类周期现象较明显. 在分子间势能、等温压缩系数和热运动效应的共同作用下, 当压力(P)为1.1倍的临界压力(Pc)时, 排列熵在0.99倍的拟临界温度(Tpc)达到最小值, P = 1.3Pc和1.5Pc时, 最小排列熵与等温压缩系数的最大值工况点保持一致, 压力继续增大, 各模拟工况密度和排列熵的波动减弱, 流体均匀性增强.Supercritical fluids (SCF) have been widely utilized in the industrial processes, such as extraction, cleaning, drying, foaming and power generation driven by primary energy. Therefore, SCF have attracted more and more attention in recent years. At supercritical state, liquid, and gas phase are not clearly distinguished, but the thermal-physical properties of fluid show an interesting characteristic, especially near the pseudo-critical temperature. Thus, it is of great significant to study the structure and density time series evolution of SCF.Due to high pressure and temperature for SCF, it can be challenging to collect experimental data of SCF. However, the advantage of molecular dynamics simulation in convenience, safty and cost over experiments. Therefore, in this paper,molecular dynamics simulation was performed to investigate the fluid structure and density series fluctuation curves at supercritical state, and the influence of parameters varitation including pressure and temperature onstructural characteristics was analyzed. In the simulation system, more than 104 atoms and simple Lennard-Jones(LJ) supercritical fluids were contained. The radial distribution function(RDF), coordination number(CN), density time series curve and permutation entropy of fluids at different pressures and temperatures were calculated. At specified pressure, the position of the first peak value of RDF gradually moves to the right with the increase of temperature, and the trend weakens with the increase of pressure. CN shows a downward trend with the increase of pressure and the CN difference at different temperatures gradually decreases. Simultaneously, the CN distribution area becomes narrow with the increase of pressure. The high/low density region calibrated by CN is stable, concentrated and large area distribution at low pressure, and the average density region is small, with the increase of pressure, the area of high/low density region is only a size of a few molecular and fluctuates sharply with time, and the area of average region is constantly expanding. At relatively low pressure, the density time series curve shows the characteristic that both the fluctuation range and quasi-period are large at pseudo-critical temperature. Simultaneously, the permutation entropy obtained from the time series curve shows three cases: (i) at low pressure (P = 1.1Pc), the minimum permutation entropy is obtained under the temperature that is lower than pseudo-critical temperature, and the system has higher orderliness; (ii) at moderate pressure (P = 1.3Pc and 1.5Pc), the state points corresponding to minimum permutation entropy is consistent with that corresponding to the maximum of isothermal compression coefficient and (iii) at high pressure (P = 2.0Pc), the permutation entropy curve fluctuates slightly and remains basically on the horizontal line. The results provide reliable support for revealing the characteristics of SCF from the microscale, and also provide useful inspiration for the practical application of SCF.
-
Keywords:
- supercritical fluids /
- molecular dynamics /
- structure characteristics /
- permutation entropy
1. 引 言
严格意义上超临界流体是指温度和压力均在临界点以上的流体[1,2], 传统教科书认为超临界态是一个均匀状态, 没有气液区分. 然而, 随着超临界流体在工业中广泛的应用, 越来越多的研究得到开展. 已有研究[3,4]发现, 可以根据定压比热容极值点的位置将超临界区分为类液和类气区, 即Widom线, 利用速度自相关函数的突变点的连线, 得到的Frenkel线是类液类气区分界线的另一个划分标准[5,6]. 此外, Banuti等[7,8]采用分子动力学和理论推导的方法进一步将超临界流体划分为类液、类气和过渡三个区域, 而分子动力学模拟的方法以其成本低、安全、便捷等优点在超临界流体模拟中得到广泛的应用.
目前关于超临界流体的分子动力学模拟主要集中在热物性和输运性质的计算以及不同势函数、截断半径选取对计算结果的影响等方面[9-12], 为后续相关研究提供可靠的设置参数. Skarmoutsos和Samios[13,14]采用分子动力学方法研究超临界水沿等温线, 系统温度(
$T$ )和临界温度(${T_{\rm{c}}}$ )的比值$T/{T_{\rm{c}}} = 1.03$ , 局部密度增强和动力学特性与系统密度的依赖关系, 此外还对该等温线上甲醇和CO2的密度增强效应进行研究, 发现两种流体均在0.7倍的临界密度(${\rho _{\rm{c}}}$ )时得到最大的增强效应; Yoshii和Okazaki[15]对Lennard-Jones (LJ) 流体沿等温线$T/{T_{\rm{c}}} = 1.03$ , 多个密度下超临界流体的结构和团簇进行研究, 认为静态结构因子在临界区域内存在较强的发散行为, 临界密度附近存在临界慢化现象[16]; Metatla等[17]根据径向分布函数、计算配位数与期望配位数之间的关系研究确定400 ℃超临界水是存在高/低密度区的非均质结构; Maddox等[18]采用分子动力学方法模拟了二维LJ流体在接近临界点时的密度不均匀特性, 指出该现象主要是“势能诱导”和“临界波动”两种作用机制共同影响的结果; Yamane等[19]不但分析了径向分布函数和静态结构因子分布曲线与温度的关系, 而且研究了近临界点附近流体汞的配位数和截断半径的依赖关系, 得到在近临界点附近流体结构分子动力学模拟需要使用大规模的系统和较大的截断半径. 除分子动力学外还有部分学者采用拉曼散射和小角度X射线散射的实验方法对超临界流体的非均匀性进行大量的研究[20-23].综上所述, 现有文献对超临界流体的结构特性系统的研究较少, 文献几乎都集中在单一的等温线上, 涉及到状态点较少, 在相图上的位置不明确, 而且没有考虑在跨越Widom线前后的结构特性的转变. 即在超临界压力下, 拟临界点附近流体的结构特性和随时间演化过程的研究仍处于空白. 因此, 针对超临界压力拟临界点附近超临界流体的结构和时变特性的研究是非常必要的. 本文利用分子动力学方法模拟研究不同压力、不同温度下超临界流体的结构特点及空间和时间演化特性, 分析压力、温度对相关物理参量的影响, 首次在分子动力学模拟中引用排列熵的概念对密度时间序列曲线进行分析, 剖析不同工况下高/低密度区形成及产生最小排列熵的主要作用机制. 研究结果为从微尺度层面揭示超临界流体的特性提供可靠支撑, 也对超临界流体的实际应用提供有益启发. 本模拟采用开源分子动力学模拟软件LAMMPS 实现, 分子位型采用Ovito软件进行可视化.
2. 物理模型及模拟细节
图1(a)为物理模型示意图. 系统沿着x, y, z三个方向均采用周期性边界条件. 模拟体系尺寸为
${L_x} = {L_y} = {L_z} = L$ , 模拟系统内充满超临界流体氩, 初始氩原子按照FCC晶格排列方式布置. 有文献[19]研究表明, 超临界流体模拟过程中需要较大的系统尺寸, 得到的计算结果具有较高的可靠性, 因此本文模拟过程中, 各模拟工况均维持系统原子数为104个, 根据不同计算工况的温度和压力, 确定各工况对应的密度, 得到模拟体系的尺寸范围为27.2283$\sigma$ —40.8040$\sigma $ .
图 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.超临界流体氩的临界点温度(Tc)为150.687 K, 压力(
${P_{\rm{c}}}$ )为4.863 MPa, 密度(${\rho _{\rm{c}}}$ )为535.6 kg/m3, 为揭示不同的超压力、拟临界温度(Tpc)附近超临界流体结构及密度时间序列曲线波动特性, 选择1.1${P_{\rm{c}}}$ —2.0${P_{\rm{c}}}$ 4个超临界压力, 如图1(b)所示, 每个压力下取0.95${T_{{\rm{pc}}}}$ —1.05${T_{{\rm{pc}}}}$ 7个工况进行模拟分析.根据超临界流体的物性参数和相关研究结果可以得到, 定压比热容(cp)的极值点连成的线即为Widom线, 也是超临界区类液类气的分界曲线, 计算工况在定压比热容的曲线上的分布如图2(a)所示, 计算工况包括了类液和类气区域, 且都集中在拟临界点附近. 相应的等温压缩系数(kT)的变化曲线如图2(b)所示, 从图中可以得到在
$P = 1.1{P_{\rm{c}}}$ 时, 等温压缩系数在拟临界温度下取得极值, 随着压力的增大, 取得极值点的位置发生右偏, 在$P = 1.3{P_{\rm{c}}}$ 和$1.5{P_{\rm{c}}}$ 时, 均在$1.01{T_{{\rm{pc}}}}$ 得到极值点, 当压力增大到$2.0{P_{\rm{c}}}$ 时, 等温压缩系数在文中选取的计算温度区间内呈现一个平缓的上升趋势, 没有出现极值点. 超临界流体物性突变是流体结构发生变化的具体表现, 也是各工况计算分析中需要重点考虑的因素.
图 2 (a) 定压比热容(cp)变化曲线; (b) 等温压缩系数(kT)变化曲线Fig. 2. (a) The curve of cp under different pressures; (b) the curve of kT under different pressures.分子动力学的基础是牛顿第二定律, 直接给出了加速度和施加外力之间的关系:
$${F_i} = m{a_i} = m\frac{{{\rm{d}}{{ r}_i}}}{{{\rm{d}}{t^2}}},$$ (1) 其中m和
$ r$ 分别为分子的质量和矢量位置, 下标i表示原子i, 右端的二阶位移导数项表示原子i的加速度${a_i}$ ,${F_i}$ 为原子i所受的总作用力.流体氩分子之间的相互作用采用Lennard-Jones (LJ) 势能模型, 表达式为
$$\phi (r) = 4\varepsilon \left[ {{{\left( {\frac{\sigma }{r}} \right)}^{12}} - {{\left( {\frac{\sigma }{r}} \right)}^6}} \right],$$ (2) 其中r为原子对间的距离, 液体氩原子之间的尺寸参数
$\sigma=0.3405\;{\rm{nm}} $ , 能量参数ε = 1.67 × 10–21J, 分子质量m = 6.69 × 10–23g. 根据相关文献研究, 在超临界模拟过程中势能截断半径为5.88$\sigma $ [24], 超过此距离的分子, 其相互作用忽略不计.采用Velocity-Verlet算法求解运动方程, 时间步长取为1 fs. 在模拟过程中对整个系统施加NVT系综, 采用Nose-Hoover控温方法. 每个算例计算15 ns, 前10 ns为充分弛豫平衡过程, 后5 ns用于统计分析各种参量, 后续的统计过程中, 忽略弛豫过程的时间, 直接从统计时刻开始记录时间为0τ.
3. 结果与讨论
3.1 径向分布函数及配位数
径向分布函数(RDF)是系统的区域密度与平均密度之比, 分子动力学中计算径向分布函数的方法为[25]
$$g({r_{\rm{c}}}) = \frac{{\rm{1}}}{{\rho 4{\text{π}}r_{\rm{c}}^2{{\text{δ}}}r}}\frac{{\displaystyle\sum\limits_{t = 1}^D {\sum\limits_{j = 1}^N {{{\Delta }}N({r_{\rm{c}}} \to {r_{\rm{c}}} +{\text{δ}}r)} } }}{{N \times D}},$$ (3) 其中
$\rho $ 为系统的密度, N为分子的总数目, D为计算的总时间(步数),${\text{δ}}r$ 为设定的距离差, 参考分子与其中心的距离${r_{\rm{c}}} \to {r_{\rm{c}}} + {\text{δ}} r$ 间的分子数目为$\Delta N$ .由图3得知, 在各计算工况下, 温度低于拟临界点温度时, 径向分布函数会出现高低不同的峰谷, 在1附近产生振荡, 表现为一种“短程有序、长程无序”的类液体现象, 和文献[26]中描述的液相趋势一致. 随着温度的升高, 波谷逐渐抬升, 第二波峰逐渐变小, 在
$T = 1.05{T_{{\rm{pc}}}}$ 时第二个峰值基本消失, 满足气相在短程区域有一个峰值, 而后单调下降到1的变化趋势, 在各模拟工况下, 随温度的升高径向分布函数表现为类液-类气的特征, 两种区域的结构存在差异. 在定压下, 随着温度的升高, 密度逐渐降低, 第一峰值出现的位置逐渐向右偏移, 随压力的增大, 偏移程度逐渐减小, 在足够大的压力下该偏移量将会消失. 四个分图对比发现, 随着压力的增大, 第一峰值逐渐减小, 峰宽度也减小, 这主要是由于超临界流体氩局部出现聚集的高密度区的结果. 氩的局部聚集必然会导致体系内部出现“空隙”, 压力的增大将这些“空隙”压缩, 减小了局部聚集氩相对密度, 因而会出现峰值和峰宽度都减小的现象, 和现有文献[27]分析保持一致. 同时也进一步证明超临界流体的高压缩性. 在$P = 1.1{P_{\rm{c}}}$ 和$1.3{P_{\rm{c}}}$ 时, 随着温度的升高峰值呈现出先下降、而后升高再下降的一种变化规律, 在$T = 1.0{T_{{\rm{pc}}}}$ 和$1.01{T_{{\rm{pc}}}}$ 时得到最大值; 在$P = 1.5{P_{\rm{c}}}$ 和$2.0{P_{\rm{c}}}$ 时, 随着温度的增加, 峰值整体呈现一个下降趋势, 在$T = 1.01{T_{{\rm{pc}}}}$ 时仍呈现微弱的增大趋势, 但是峰值最大值仍出现在最低温度工况下, 这主要由于在拟临界点附近超临界流体存在很大的压缩性, 促使流体聚集产生高密度区, 但是随着压力的增加, 流体的压缩性迅速降低, 温度升高带来的分子热运动逐渐加强, 使得流体的聚集能力逐渐减弱.
图 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.配位数(CN)的变化趋势和密度呈线性关系, 也是描述流体微观结构的重要物理参数, 反映了距离中心分子为rx的球体区域内分子的个数, 描述了分子排列的紧密程度, 配位数越大, 粒子排列越紧密, 不同区间范围内的配位数的定义计算式为[28]
$${N_{\rm{c}}} = \int\nolimits_0^{{r_x}} {g({r_{\rm{c}}})} 4{\text{π}}\rho r_{\rm{c}}^2{\rm{d}}{r_{\rm{c}}},$$ (4) 计算超临界流体氩的配位数能清楚反映出流体的内部结构. 具体模拟结果如图4所示.
图 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.有研究表明配位数不仅与结构、密度有关, 而且是温度的函数[29]. 在相同的压力下, 随着温度的升高, 流体的密度减小, 此时分子间的距离增大, 分子间引力对配位数的影响越来越大, 但是温度的升高, 促使分子的热运动加剧, 在引力和热运动的共同作用下, 发现随温度增加配位数逐渐减小. 此外随着压力的增大, 确定的温度区间内, 配位数的波动范围逐渐缩小, 各工况的差异性逐渐减弱. 从图5可以观察到,
$P = 1.1{P_{\rm{c}}}$ 时, 随温度的增加, 配位数曲线的斜率在16.22—5.90区间变化, 随着压力的升高;$P = 2.0{P_{\rm{c}}}$ 时, 曲线斜率在14.27—9.25区间变化. 压力越接近临界点时, 随着温度的变化, 密度波动范围越大, 从而引起配位数的波动区间较大, 而在远离临界点的压力下, 密度随温度的变化范围较小, 相应的配位数曲线的斜率范围较窄. 在$T/{T_{{\rm{pc}}}} < 1.0$ 的类液区, 压力增大, 拟临界点温度升高, 密度减小, 配位数斜率呈现减小的变化趋势; 在$T/{T_{{\rm{pc}}}} \geqslant 1.0$ 的类气区表现为相反的变化趋势. 可以得到密度是影响流体配位数的主要参数, 温度和压力也是通过调整密度的大小间接影响配位数分布.根据配位数的计算值(
${N_{\rm{c}}}$ )和期望值(${N_{\rm{e}}}$ )对密度不均匀性进行判断, 如果直接比较两个数字的大小, 这样的标准过于严格, 局部结构的微小波动, 即使只有一个分子穿过设置的半径边界, 也会导致密度的分类从平均转变成高或低密度. 根据文献[30]提到的稍微宽松且能准确判断的标准, 即选取一个小量的值$\delta $ , 用${N_{\rm{e}}} \pm \delta $ 作为参数度量平均密度的变化, 具体划分原则如下:$$ \begin{split} & {\rm{high\; if}}\;{N_{\rm{c}}} > {N_{\rm{e}}} + \delta,\\[2mm] &\quad {\rm{average \;if }}\;{N_{\rm{e}}} - \delta \leqslant {N_{\rm{c}}} \leqslant {N_{\rm{e}}} + \delta,\\[2mm] & {\rm{low\; if}}\;{N_c} < {N_{\rm{e}}} - \delta . \end{split} $$ (5) 计算中允许的波动量为30%
${N_{\rm{e}}}$ , 则$\delta $ 应满足2$\delta $ +1=0.3${N_{\rm{e}}}$ , 进一步得到不同参数下$\delta $ 的具体值, 利用划分原则可以判断密度分布趋势, 称该方法为“30%方法”. 采用该方法, 根据配位数得到在$P = 1.1{P_{\rm{c}}}$ 和$2.0{P_{\rm{c}}}$ ,$T = {T_{{\rm{pc}}}}$ 流体在xy平面内高/低及平均密度区分布随时间的演化如图6所示.
图 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.图6(a)表示在0—5000τ的时间范围内
$P = 1.1{P_{\rm{c}}}$ ,$T = {T_{{\rm{pc}}}}$ 时的演化过程, 由图可知, 该计算工况产生的均值区面积较小, 低密度区的位置主要集中在模拟系统的中部, 随着时间的演化, 密度在不停地波动, 低密度区呈现分散-聚合-分散的演化规律, 产生的高密度区的面积相对较大且位置集中, 演化过程中波动微弱. 由于在超临界区表面张力的消失, 不存在亚临界工况下的弯曲界面. 形成该现象的原因主要是由于当密度低时, 分子间引力起主导作用, 在温度的影响下分子处于不断的相互碰撞中, 能量的增加, 导致高密度区的形成. 此外, 该工况拟临界点位置出现等温压缩系数的极值点, 形成较高的密度区所付出的代价变小, 超临界流体的关联长度在拟临界处也存在极大值, 因此会出现大面积的, 相对稳定的高密度区. 在确定尺寸和分子数的系统中, 高密度区的形成必然会引起低密度区的产生. 随着高密度区的形成, 分子间的距离缩短, 增大的斥力将部分分子从高密度区向低密度区推, 当分子间距离较大时, 引力起主导作用, 会使得分子再次聚集成高密度区, 但此时作用势效应相对等温压缩系数效应较小, 因此各区域所处位置稳定, 波动微弱. 随着压力的升高, 流体的等温压缩系数仍存在极值点, 但是数值较小, 流体可压缩性减弱, 较难形成高密度区, 仅在分子间短程势作用下, 形成由几个分子组成的且较为分散的高密度区, 形成相对较大的平均值区, 如图6(b)所示. 随时间的演化, 各区域不停波动, 存在明显的嵌套现象, 系统中局部出现类似“花斑”的现象, 这些花斑若隐若现、此起彼伏、互相嵌套的性质和相关教课书[31]中提出的结论保持一致.从图7可以直观地观察各工况不同密度区所占比例随时间的变化. 不同的压力条件下, 类液和类气区的占比一直处于一个波动的状态, 仅观察曲线, 发现波动过程并没有实际的规律可循, 整体表现为一个混乱的动态变化过程. 从图7(a)和图7(b)分图中发现随着压力增大类液和类气区的占比整体呈现下降的趋势, 进一步说明随压力增大, 形成高密度区较为困难, 计算压力距临界压力越远, 流体的均匀性越强. 在
$P = 1.1{P_{\rm{c}}}$ ,$T = {T_{{\rm{pc}}}}$ 的工况时, 高密度区占比约为60%, 低密度区的占比约为35%, 此时均匀区域的占比最小, 系统的不均匀性最强.
图 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.3.2 静态结构因子
结构因子主要表征材料对射线的散射能量, 反映材料结构的平均信息, 可以进一步应用到流体中, 观察流体的结构特性. 而静态结构因子和径向分布函数互为傅里叶变换[19,32], 计算式为
$$S(k) = 1 + n\int {\left[ {g({r_{\rm{c}}}) - 1} \right]} {\kern 1pt} {{\rm{e}}^{{\rm{i}}k \cdot {r_{\rm{c}}}}}{\rm{d}}{r_{\rm{c}}},$$ (6) 其中n = N/V, V为系统体积, 由于结构因子是倒易空间, 变量k与距中心分子的距离
${r_{\rm{c}}}$ 成反比.由图8可知, 各模拟工况下均在较小的k值范围内存在发散行为, k < 0.5
${\sigma ^{ - 1}}$ , 随着压力的增加, 这种发散行为逐渐减弱. 在$P = 1.1{P_{\rm{c}}}$ 时, 拟临界点处的发散行为最为强烈, 流体表现为较强的小角度散射, 曲线在一定范围内存在微小的波动, 而随着压力的升高, 这种波动现象消失. 在$P = 1.3{P_{\rm{c}}}$ 和$1.5{P_{\rm{c}}}$ 时, 发散最强烈的行为出现在偏离拟临界点的$T = 1.01{T_{{\rm{pc}}}}$ 工况, 与等温压缩系数的最大值点保持一致. 在图8(a)—图8(c)分图中得到随着温度偏离拟临界点的距离增大, 发散性减小, 即各工况均在$T = 0.95{T_{{\rm{pc}}}}$ 和$1.05{T_{{\rm{pc}}}}$ 时得到静态结构因子的最小值. 在图8(d)分中, 由于等温压缩系数在一个较小的水平缓慢增加, 因此静态结构因子的发散特性也呈现微弱增加趋势. 压力较低时, 静态结构因子的发散在拟临界点温度达到最大值, 但随着压力的增加, 发散的极值点工况与等温压缩系数的极值点工况相符合.
图 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.3.3 密度时间序列曲线和排列熵
密度随时间的变化是分子动力学模拟中比较容易得到的数据, 该曲线给出初步粗略的波动信息. 从图9(a)可以看出, 压力
$P = 1.1{P_{\rm{c}}}$ , 密度时间序列曲线呈现出包含类周期特征的大幅波动, 流体结构具有局部有序特征, 但这种波动存在不规则运动的相互叠加, 大波套小波的现象, 并不具有严格的周期性. 随着压力的增加, 波动的幅度和类周期性均减弱, 随机性增强, 此时流体结构具有较强的无序性.
图 9 (a) 密度时间序列曲线 ( T = Tpc); (b) 排列熵Fig. 9. (a) Time series of density for T = Tpc; (b) permutation entropy.为具体描述密度时间序列的复杂程度, 引入排列熵的分析方法, 对一个确定长度的时间序列, 可以根据自相关函数法确定一个延迟时间
${\tau _{\rm{D}}}$ , 根据排列熵嵌入维度的要求, 文中选取各工况嵌入维数m = 5. 排列熵作为衡量时间序列曲线复杂度的指标, 越规则的时间序列, 对应的排列熵越小, 越复杂的时间序列, 对应的排列熵越大.根据文献[33]提出的计算方法, 对各工况的密度时间序列曲线进行排列熵的计算, 从图9(b)中可以观察到在
$P = 1.1{P_{\rm{c}}}$ 时, 在$T = 0.99{T_{{\rm{pc}}}}$ 和$1.0{T_{{\rm{pc}}}}$ 时排列熵的值较小, 对应较规则的时间序列. 随着压力的增大,$P = 1.3{P_{\rm{c}}}$ 和$1.5{P_{\rm{c}}}$ 时, 均在$T = 1.01{T_{{\rm{pc}}}}$ 时得到最小的排列熵, 随着压力的进一步增加, 排列熵的值呈微弱增大趋势, 且各工况的值在一个水平线附近波动. 产生这种现象的原因主要是分子间作用势、等温压缩系数和分子热运动效应共同作用的结果.$P = 1.1{P_{\rm{c}}}$ 时, 各计算工况的温度低, 分子热运动效应较弱, 等温压缩系数在拟临界点处达到极大值, 流体容易被压缩形成高密度区, 在这两种效应的共同作用下, 排列熵的最小值出现在$T = 0.99{T_{{\rm{pc}}}}$ 的工况. 随着压力的升高($P = 1.3{P_{\rm{c}}}$ 和$1.5{P_{\rm{c}}}$ ), 拟临界温度逐渐增大, 分子的热运动加剧, 导致分子间的作用势效应减小, 此时等温压缩系数效应在密度波动中占据主导地位, 最小排列熵和等温压缩系数极大值点工况一致. 当$P = 2.0{P_{\rm{c}}}$ 时, 各工况得到的排列熵变化趋势和等温压缩系数保持一致.4. 结 论
采用分子动力学方法模拟不同超临界压力, 拟临界温度附近流体的结构特征, 对流体的径向分布函数、配位数、静态结构因子、密度时间序列曲线及排列熵展开研究, 分析压力和温度的变化对各参量的影响, 得到如下结论:
1) 径向分布函数在类液区存在“短程有序、长程无序”的现象, 随着温度的升高, 这种有序的现象逐渐消失, 在类气区仅在短程区域存在一个峰值, 而后单调下降到; 定压下, 随着温度的升高, 第一峰值出现的位置逐渐右移, 但这种变化趋势随着压力偏离临界点距离的增大而减弱; 此外首次提出, 在
$P = 1.1{P_{\rm{c}}}$ 和$1.3{P_{\rm{c}}}$ 时, 峰值的最大值均在等温压缩极大值点工况出现, 但是随着压力的增加, 该现象不再明显.2) 定压时, 随着温度的升高, 配位数逐渐减小, 压力的增加导致不同温度下流体密度的差值减小, 因此引起配位数的波动范围缩小, 分布区域变窄; 采用配位数标定高/低密度区的标准, 在
$P = 1.1{P_{\rm{c}}}$ ,$T = {T_{{\rm{pc}}}}$ 时可以得到大面积、分布集中且波动较小的高/低密度区, 此时均值区域占比较小; 随着压力升高, 均值区域逐渐增加, 高/低密度区逐渐减小, 仅有几个分子大小, 且相互嵌套, 并随时间有较明显的波动, 说明随着压力的增大, 流体的均匀性逐渐增强.3) 静态结构因子是通过散射效应观察流体的内部结构, 各模拟工况均在较小的k值内存在发散行为, 随着压力的增加, 静态结构因子曲线的发散值呈迅速下降的趋势, 研究发现定压时静态结构因子最大值和等温压缩系数极值点工况保持一致.
4) 密度时间序列曲线可以初步得到中心切片密度随压力的增加, 波动幅度逐渐减弱的规律, 和前文压力增大流体均匀性增强的结论符合. 排列熵的变化规律可以分为三类: 一是在
$P = 1.1{P_{\rm{c}}}$ 的低压下, 排列熵在小于拟临界温度的位置得到最大值; 二是在$P = 1.3{P_{\rm{c}}}$ 和$1.5{P_{\rm{c}}}$ 的中压下, 排列熵的最小值点和等温压缩系数极值点一致; 三是在P =$2.0{P_{\rm{c}}}$ 的高压下, 排列熵的值在一个水平线附近波动, 整体变化趋势较为平缓, 流体均匀性增强. 探讨微尺度下超临界流体的结构特征有助于全面了解超临界流体特性, 为超临界流体的应用提供有力支撑.[1] Pena-Pereira F, Tobiszewski M 2017 Elsevier 155
[2] Carlès P 2010 J. Supercrit. Fluids 53 2
Google Scholar
[3] Raju M, BanutiD T, MaP C, Ihme M 2017 Sci. Rep.-UK 7 3027
Google Scholar
[4] Artemenko S, Krijgsman P, Mazur V 2017 J. Mol. Liq. 238 122
Google Scholar
[5] Brazhkin V V, Fomin Y D, Lyapin A G, Lyapin V V, Ryzhov V N, TsiokE N, Trachenko Kostya 2013 Phys. Rev. Lett. 111 145901
Google Scholar
[6] Fomin Y D, Ryzhov V N, Tsiok E N, Brazhkin V V 2015 Sci.Rep.-UK 5 14234
Google Scholar
[7] Banuti D T, Raju M, Ihme M 2017 Cent. Turbul. Res. Annu. Res. Briefs 165
[8] Banuti D T 2015 J.Supercrit.Fluids. 98 12
Google Scholar
[9] Raman A S, Li H, Chiew Y C 2018 J. Chem. Phys. 148 014502
Google Scholar
[10] Nichele J, Abreu C R A, Alves L S B, Jr I B 2018 J.Supercrit.Fluids. 135 225
Google Scholar
[11] Nichele J, de Oliveira A B, Alves L S B, Borges I 2017 J. Mol.Liq. 237 65
Google Scholar
[12] Egorov S A 2002 Chem. Phys. Lett. 354 140
Google Scholar
[13] Skarmoutsos I, Samios J 2007 J. Chem. Phys. 126 044503
Google Scholar
[14] Skarmoutsos I, Samios J 2006 J. Phys. Chem. B 110 21931
Google Scholar
[15] Yoshii N, Okazaki S 1998 Fluid PhaseEquilib. 144 225
[16] Yoshii N, Okazaki S 1997 J. Chem. Phys. 107 2020
Google Scholar
[17] Metatla N, Lafond F, Jay-Gerin J P, Soldera A 2016 Rsc. Adv. 6 30484
Google Scholar
[18] Maddox M W, Goodyear G, Tucker S C 2000 J. Phys. Chem. B 104 6248
Google Scholar
[19] Yamane A, Shimojo F, Hoshino K 2006 J. Phys. Soc.Jpn. 75 124602
Google Scholar
[20] Nishikawa K, Arai A A, Morita T 2004 J.Supercrit.Fluids. 30 249
Google Scholar
[21] Nishikawa K, Ochiai H, Saitow K, Morita T 2003 Chem. Phys. 286 421
Google Scholar
[22] Cabaço M I, Besnard M, Tassaing T, Danten Y 2004 Pure Appl. Chem. 76 141
Google Scholar
[23] Arai A A, Morita T, Nishikawa K 2007 Fluid Phase Equilibria. 252 114
Google Scholar
[24] Ghosh K, Krishnamurthy C V 2018 Phys. Rev. E 97 012131
[25] 陈正隆, 徐为人, 汤立达 2007 分子模拟的理论与实践 (北京: 化学工业出版社) 第110−112页
Chen Z L, Xu W R, Tang L D 2007 (Beijing: Chemical Industry Press) pp110−112 (in Chinese)
[26] Bolmatov D, Brazhkin V V, Fomin Y D, Ryzhov V N, Trachenko K 2013 J. Chem. Phys. 139 234501
Google Scholar
[27] 吴方棣, 郑辉东, 刘俊劭, 郑细鸣 2014 辽宁石油化工大学学报 34 8
Google Scholar
Wu F K, Zheng H D, Liu J X, Zheng X M 2014 Journal of Liaoning Shihua University 34 8
Google Scholar
[28] Skarmoutsos I, Guardia E, Samios J 2017 J.Supercrit.Fluids. 130 156
Google Scholar
[29] 计伟荣 1993 浙江工业大学学报 59 1
Ji W R 1993 Journal of Zhejiang University of Technology 59 1
[30] Martinez H L, Ravi R, Tucker S C 1996 J. Chem. Phys. 104 1067
Google Scholar
[31] 于渌, 郝柏林, 陈晓松 2016 边缘奇迹:相变和临界现象 (北京: 科学出版社) 第81−93页
Yu L, Hao B L, Chen X S 2016 (Beijing: Science Press) pp81−93 (in Chinese)
[32] March N H, Tosi M P 2002 (Singapore: World Scientific) pp75−80
[33] Bandt C, Pompe B 2002 Phys. Rev. Lett. 88 174102
Google Scholar
-
图 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.
图 3 径向分布函数 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc
Fig. 3. Radial distribution function: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.
图 4 配位数 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc
Fig. 4. Coordination number: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.
图 6 流体在xy平面内高/低密度区分布随时间的演化 (a) P = 1.1Pc, T = Tpc; (b) P = 2.0Pc, T = Tpc
Fig. 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.
图 8 静态结构因子 (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc
Fig. 8. Static structure fator: (a) P = 1.1Pc; (b) P = 1.3Pc; (c) P = 1.5Pc; (d) P = 2.0Pc.
-
[1] Pena-Pereira F, Tobiszewski M 2017 Elsevier 155
[2] Carlès P 2010 J. Supercrit. Fluids 53 2
Google Scholar
[3] Raju M, BanutiD T, MaP C, Ihme M 2017 Sci. Rep.-UK 7 3027
Google Scholar
[4] Artemenko S, Krijgsman P, Mazur V 2017 J. Mol. Liq. 238 122
Google Scholar
[5] Brazhkin V V, Fomin Y D, Lyapin A G, Lyapin V V, Ryzhov V N, TsiokE N, Trachenko Kostya 2013 Phys. Rev. Lett. 111 145901
Google Scholar
[6] Fomin Y D, Ryzhov V N, Tsiok E N, Brazhkin V V 2015 Sci.Rep.-UK 5 14234
Google Scholar
[7] Banuti D T, Raju M, Ihme M 2017 Cent. Turbul. Res. Annu. Res. Briefs 165
[8] Banuti D T 2015 J.Supercrit.Fluids. 98 12
Google Scholar
[9] Raman A S, Li H, Chiew Y C 2018 J. Chem. Phys. 148 014502
Google Scholar
[10] Nichele J, Abreu C R A, Alves L S B, Jr I B 2018 J.Supercrit.Fluids. 135 225
Google Scholar
[11] Nichele J, de Oliveira A B, Alves L S B, Borges I 2017 J. Mol.Liq. 237 65
Google Scholar
[12] Egorov S A 2002 Chem. Phys. Lett. 354 140
Google Scholar
[13] Skarmoutsos I, Samios J 2007 J. Chem. Phys. 126 044503
Google Scholar
[14] Skarmoutsos I, Samios J 2006 J. Phys. Chem. B 110 21931
Google Scholar
[15] Yoshii N, Okazaki S 1998 Fluid PhaseEquilib. 144 225
[16] Yoshii N, Okazaki S 1997 J. Chem. Phys. 107 2020
Google Scholar
[17] Metatla N, Lafond F, Jay-Gerin J P, Soldera A 2016 Rsc. Adv. 6 30484
Google Scholar
[18] Maddox M W, Goodyear G, Tucker S C 2000 J. Phys. Chem. B 104 6248
Google Scholar
[19] Yamane A, Shimojo F, Hoshino K 2006 J. Phys. Soc.Jpn. 75 124602
Google Scholar
[20] Nishikawa K, Arai A A, Morita T 2004 J.Supercrit.Fluids. 30 249
Google Scholar
[21] Nishikawa K, Ochiai H, Saitow K, Morita T 2003 Chem. Phys. 286 421
Google Scholar
[22] Cabaço M I, Besnard M, Tassaing T, Danten Y 2004 Pure Appl. Chem. 76 141
Google Scholar
[23] Arai A A, Morita T, Nishikawa K 2007 Fluid Phase Equilibria. 252 114
Google Scholar
[24] Ghosh K, Krishnamurthy C V 2018 Phys. Rev. E 97 012131
[25] 陈正隆, 徐为人, 汤立达 2007 分子模拟的理论与实践 (北京: 化学工业出版社) 第110−112页
Chen Z L, Xu W R, Tang L D 2007 (Beijing: Chemical Industry Press) pp110−112 (in Chinese)
[26] Bolmatov D, Brazhkin V V, Fomin Y D, Ryzhov V N, Trachenko K 2013 J. Chem. Phys. 139 234501
Google Scholar
[27] 吴方棣, 郑辉东, 刘俊劭, 郑细鸣 2014 辽宁石油化工大学学报 34 8
Google Scholar
Wu F K, Zheng H D, Liu J X, Zheng X M 2014 Journal of Liaoning Shihua University 34 8
Google Scholar
[28] Skarmoutsos I, Guardia E, Samios J 2017 J.Supercrit.Fluids. 130 156
Google Scholar
[29] 计伟荣 1993 浙江工业大学学报 59 1
Ji W R 1993 Journal of Zhejiang University of Technology 59 1
[30] Martinez H L, Ravi R, Tucker S C 1996 J. Chem. Phys. 104 1067
Google Scholar
[31] 于渌, 郝柏林, 陈晓松 2016 边缘奇迹:相变和临界现象 (北京: 科学出版社) 第81−93页
Yu L, Hao B L, Chen X S 2016 (Beijing: Science Press) pp81−93 (in Chinese)
[32] March N H, Tosi M P 2002 (Singapore: World Scientific) pp75−80
[33] Bandt C, Pompe B 2002 Phys. Rev. Lett. 88 174102
Google Scholar
下载:
下载:
计量
- 文章访问数: 12061
- PDF下载量: 187
- 被引次数: 0
