超臨界Lennard-Jones流體結(jié)構(gòu)特性分子動力學研究*

王艷 徐進良2)? 李文 劉歡

1) (華北電力大學, 低品位能源多相流與傳熱北京市重點實驗室, 北京 102206)

2) (華北電力大學, 電站能量傳遞轉(zhuǎn)化與系統(tǒng)教育部重點實驗室, 北京 102206)

研究超臨界流體在不同壓力和溫度的結(jié)構(gòu)特征有助于深刻理解并有效利用超臨界流體.本文采用分子動力學方法模擬超臨界壓力、擬臨界溫度附近流體的結(jié)構(gòu)及密度波動曲線的排列熵, 分析狀態(tài)參數(shù)變化的影響.結(jié)果表明, 定壓下, 徑向分布函數(shù)隨溫度升高, 第一峰值位置逐漸向右移動, 但右移幅度隨著壓力偏離臨界點距離的增大而減弱, 近臨界壓力時, 出現(xiàn)峰值最高點的工況和等溫壓縮系數(shù)的極值點位置一致, 壓力增大, 該現(xiàn)象消失.低壓力擬臨界點時易出現(xiàn)面積大、相對集中且分布穩(wěn)定的高/低密度區(qū), 無明顯嵌套現(xiàn)象.靜態(tài)結(jié)構(gòu)因子存在一定發(fā)散行為, 發(fā)散的最大值和等溫壓縮系數(shù)極值點所處工況符合.低壓力時密度時間序列的波動幅度最大, 類周期現(xiàn)象較明顯.在分子間勢能、等溫壓縮系數(shù)和熱運動效應(yīng)的共同作用下, 當壓力(P)為 1.1 倍的臨界壓力 (Pc)時, 排列熵在 0.99 倍的擬臨界溫度 (Tpc)達到最小值, P = 1.3Pc和 1.5Pc時, 最小排列熵與等溫壓縮系數(shù)的最大值工況點保持一致, 壓力繼續(xù)增大, 各模擬工況密度和排列熵的波動減弱,流體均勻性增強.

1 引 言

嚴格意義上超臨界流體是指溫度和壓力均在臨界點以上的流體[1,2], 傳統(tǒng)教科書認為超臨界態(tài)是一個均勻狀態(tài), 沒有氣液區(qū)分.然而, 隨著超臨界流體在工業(yè)中廣泛的應(yīng)用, 越來越多的研究得到開展.已有研究[3,4]發(fā)現(xiàn), 可以根據(jù)定壓比熱容極值點的位置將超臨界區(qū)分為類液和類氣區(qū), 即Widom線, 利用速度自相關(guān)函數(shù)的突變點的連線,得到的Frenkel線是類液類氣區(qū)分界線的另一個劃分標準[5,6].此外, Banuti等[7,8]采用分子動力學和理論推導的方法進一步將超臨界流體劃分為類液、類氣和過渡三個區(qū)域, 而分子動力學模擬的方法以其成本低、安全、便捷等優(yōu)點在超臨界流體模擬中得到廣泛的應(yīng)用.

目前關(guān)于超臨界流體的分子動力學模擬主要集中在熱物性和輸運性質(zhì)的計算以及不同勢函數(shù)、截斷半徑選取對計算結(jié)果的影響等方面[9?12], 為后續(xù)相關(guān)研究提供可靠的設(shè)置參數(shù).Skarmoutsos和Samios[13,14]采用分子動力學方法研究超臨界水沿等溫線, 系統(tǒng)溫度( T )和臨界溫度( Tc)的比值T/Tc=1.03, 局部密度增強和動力學特性與系統(tǒng)密度的依賴關(guān)系, 此外還對該等溫線上甲醇和CO2的密度增強效應(yīng)進行研究, 發(fā)現(xiàn)兩種流體均在0.7倍的臨界密度( ρc)時得到最大的增強效應(yīng);Yoshii和 Okazaki[15]對 Lennard-Jones (LJ) 流體沿等溫線 T /Tc=1.03 , 多個密度下超臨界流體的結(jié)構(gòu)和團簇進行研究, 認為靜態(tài)結(jié)構(gòu)因子在臨界區(qū)域內(nèi)存在較強的發(fā)散行為, 臨界密度附近存在臨界慢化現(xiàn)象[16]; Metatla等[17]根據(jù)徑向分布函數(shù)、計算配位數(shù)與期望配位數(shù)之間的關(guān)系研究確定400 ℃超臨界水是存在高/低密度區(qū)的非均質(zhì)結(jié)構(gòu);Maddox等[18]采用分子動力學方法模擬了二維LJ流體在接近臨界點時的密度不均勻特性, 指出該現(xiàn)象主要是“勢能誘導”和“臨界波動”兩種作用機制共同影響的結(jié)果; Yamane等[19]不但分析了徑向分布函數(shù)和靜態(tài)結(jié)構(gòu)因子分布曲線與溫度的關(guān)系,而且研究了近臨界點附近流體汞的配位數(shù)和截斷半徑的依賴關(guān)系, 得到在近臨界點附近流體結(jié)構(gòu)分子動力學模擬需要使用大規(guī)模的系統(tǒng)和較大的截斷半徑.除分子動力學外還有部分學者采用拉曼散射和小角度X射線散射的實驗方法對超臨界流體的非均勻性進行大量的研究[20?23].

綜上所述, 現(xiàn)有文獻對超臨界流體的結(jié)構(gòu)特性系統(tǒng)的研究較少, 文獻幾乎都集中在單一的等溫線上, 涉及到狀態(tài)點較少, 在相圖上的位置不明確,而且沒有考慮在跨越Widom線前后的結(jié)構(gòu)特性的轉(zhuǎn)變.即在超臨界壓力下, 擬臨界點附近流體的結(jié)構(gòu)特性和隨時間演化過程的研究仍處于空白.因此, 針對超臨界壓力擬臨界點附近超臨界流體的結(jié)構(gòu)和時變特性的研究是非常必要的.本文利用分子動力學方法模擬研究不同壓力、不同溫度下超臨界流體的結(jié)構(gòu)特點及空間和時間演化特性, 分析壓力、溫度對相關(guān)物理參量的影響, 首次在分子動力學模擬中引用排列熵的概念對密度時間序列曲線進行分析, 剖析不同工況下高/低密度區(qū)形成及產(chǎn)生最小排列熵的主要作用機制.研究結(jié)果為從微尺度層面揭示超臨界流體的特性提供可靠支撐, 也對超臨界流體的實際應(yīng)用提供有益啟發(fā).本模擬采用開源分子動力學模擬軟件LAMMPS 實現(xiàn), 分子位型采用Ovito軟件進行可視化.

2 物理模型及模擬細節(jié)

圖1(a)為物理模型示意圖.系統(tǒng)沿著 x, y,z三個方向均采用周期性邊界條件.模擬體系尺寸為 Lx=Ly=Lz=L , 模擬系統(tǒng)內(nèi)充滿超臨界流體氬, 初始氬原子按照FCC晶格排列方式布置.有文獻[19]研究表明, 超臨界流體模擬過程中需要較大的系統(tǒng)尺寸, 得到的計算結(jié)果具有較高的可靠性, 因此本文模擬過程中, 各模擬工況均維持系統(tǒng)原子數(shù)為104個, 根據(jù)不同計算工況的溫度和壓力,確定各工況對應(yīng)的密度, 得到模擬體系的尺寸范圍為27.2283 σ —40.8040 σ.

超臨界流體氬的臨界點溫度(Tc)為150.687 K,壓力 ( Pc)為4.863 MPa, 密度 ( ρc)為 535.6 kg/m3,為揭示不同的超壓力、擬臨界溫度(Tpc)附近超臨界流體結(jié)構(gòu)及密度時間序列曲線波動特性, 選擇1.1 Pc—2.0 Pc4 個超臨界壓力, 如圖1(b)所示, 每個壓力下取0.95 Tpc—1.05 Tpc7個工況進行模擬分析.

圖1 (a) 物理模型圖; (b) 模擬狀態(tài)點在相圖中的分布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) 等溫壓縮系數(shù) (kT)變化曲線Fig.2.(a) The curve of cp under different pressures; (b) the curve of kT under different pressures.

根據(jù)超臨界流體的物性參數(shù)和相關(guān)研究結(jié)果可以得到, 定壓比熱容(cp)的極值點連成的線即為Widom線, 也是超臨界區(qū)類液類氣的分界曲線,計算工況在定壓比熱容的曲線上的分布如圖2(a)所示, 計算工況包括了類液和類氣區(qū)域, 且都集中在擬臨界點附近.相應(yīng)的等溫壓縮系數(shù)(kT)的變化曲線如圖2(b)所示, 從圖中可以得到在P=1.1Pc時, 等溫壓縮系數(shù)在擬臨界溫度下取得極值, 隨著壓力的增大, 取得極值點的位置發(fā)生右偏, 在P=1.3Pc和 1.5Pc時, 均在 1.01Tpc得到極值點, 當壓力增大到 2.0Pc時, 等溫壓縮系數(shù)在文中選取的計算溫度區(qū)間內(nèi)呈現(xiàn)一個平緩的上升趨勢, 沒有出現(xiàn)極值點.超臨界流體物性突變是流體結(jié)構(gòu)發(fā)生變化的具體表現(xiàn), 也是各工況計算分析中需要重點考慮的因素.

分子動力學的基礎(chǔ)是牛頓第二定律, 直接給出了加速度和施加外力之間的關(guān)系:

其中m和 r 分別為分子的質(zhì)量和矢量位置, 下標i表示原子i, 右端的二階位移導數(shù)項表示原子i的加速度 ai, Fi為原子i所受的總作用力.

流體氬分子之間的相互作用采用Lennard-Jones (LJ) 勢能模型, 表達式為

其中r為原子對間的距離, 液體氬原子之間的尺寸參數(shù) σ =0.3405 nm , 能量參數(shù) ε = 1.67 × 10–21J,分子質(zhì)量 m = 6.69 × 10–23g.根據(jù)相關(guān)文獻研究,在超臨界模擬過程中勢能截斷半徑為5.88 σ[24], 超過此距離的分子, 其相互作用忽略不計.

采用Velocity-Verlet算法求解運動方程, 時間步長取為1 fs.在模擬過程中對整個系統(tǒng)施加NVT 系綜, 采用 Nose-Hoover控溫方法.每個算例計算 15 ns, 前 10 ns為充分弛豫平衡過程, 后5 ns用于統(tǒng)計分析各種參量, 后續(xù)的統(tǒng)計過程中,忽略弛豫過程的時間, 直接從統(tǒng)計時刻開始記錄時間為 0τ.

3 結(jié)果與討論

3.1 徑向分布函數(shù)及配位數(shù)

徑向分布函數(shù)(RDF)是系統(tǒng)的區(qū)域密度與平均密度之比, 分子動力學中計算徑向分布函數(shù)的方法為[25]

其中 ρ 為系統(tǒng)的密度, N為分子的總數(shù)目, D為計算的總時間 (步數(shù)), δr 為設(shè)定的距離差, 參考分子與其中心的距離 rc→ rc+δr 間的分子數(shù)目為 ? N.

圖3 徑向分布函數(shù) (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得知, 在各計算工況下, 溫度低于擬臨界點溫度時, 徑向分布函數(shù)會出現(xiàn)高低不同的峰谷, 在1附近產(chǎn)生振蕩, 表現(xiàn)為一種“短程有序、長程無序”的類液體現(xiàn)象, 和文獻[26]中描述的液相趨勢一致.隨著溫度的升高, 波谷逐漸抬升, 第二波峰逐漸變小, 在 T =1.05Tpc時第二個峰值基本消失, 滿足氣相在短程區(qū)域有一個峰值, 而后單調(diào)下降到1的變化趨勢, 在各模擬工況下, 隨溫度的升高徑向分布函數(shù)表現(xiàn)為類液-類氣的特征, 兩種區(qū)域的結(jié)構(gòu)存在差異.在定壓下, 隨著溫度的升高,密度逐漸降低, 第一峰值出現(xiàn)的位置逐漸向右偏移, 隨壓力的增大, 偏移程度逐漸減小, 在足夠大的壓力下該偏移量將會消失.四個分圖對比發(fā)現(xiàn),隨著壓力的增大, 第一峰值逐漸減小, 峰寬度也減小, 這主要是由于超臨界流體氬局部出現(xiàn)聚集的高密度區(qū)的結(jié)果.氬的局部聚集必然會導致體系內(nèi)部出現(xiàn)“空隙”, 壓力的增大將這些“空隙”壓縮, 減小了局部聚集氬相對密度, 因而會出現(xiàn)峰值和峰寬度都減小的現(xiàn)象, 和現(xiàn)有文獻[27]分析保持一致.同時也進一步證明超臨界流體的高壓縮性.在P=1.1Pc和 1.3Pc時, 隨著溫度的升高峰值呈現(xiàn)出先下降、而后升高再下降的一種變化規(guī)律, 在T=1.0Tpc和 1.01Tpc時得到最大值; 在 P =1.5Pc和2.0Pc時, 隨著溫度的增加, 峰值整體呈現(xiàn)一個下降趨勢, 在 T =1.01Tpc時仍呈現(xiàn)微弱的增大趨勢, 但是峰值最大值仍出現(xiàn)在最低溫度工況下, 這主要由于在擬臨界點附近超臨界流體存在很大的壓縮性,促使流體聚集產(chǎn)生高密度區(qū), 但是隨著壓力的增加, 流體的壓縮性迅速降低, 溫度升高帶來的分子熱運動逐漸加強, 使得流體的聚集能力逐漸減弱.

配位數(shù)(CN)的變化趨勢和密度呈線性關(guān)系,也是描述流體微觀結(jié)構(gòu)的重要物理參數(shù), 反映了距離中心分子為rx的球體區(qū)域內(nèi)分子的個數(shù), 描述了分子排列的緊密程度, 配位數(shù)越大, 粒子排列越緊密, 不同區(qū)間范圍內(nèi)的配位數(shù)的定義計算式為[28]

計算超臨界流體氬的配位數(shù)能清楚反映出流體的內(nèi)部結(jié)構(gòu).具體模擬結(jié)果如圖4所示.

有研究表明配位數(shù)不僅與結(jié)構(gòu)、密度有關(guān), 而且是溫度的函數(shù)[29].在相同的壓力下, 隨著溫度的升高, 流體的密度減小, 此時分子間的距離增大,分子間引力對配位數(shù)的影響越來越大, 但是溫度的升高, 促使分子的熱運動加劇, 在引力和熱運動的共同作用下, 發(fā)現(xiàn)隨溫度增加配位數(shù)逐漸減小.此外隨著壓力的增大, 確定的溫度區(qū)間內(nèi), 配位數(shù)的波動范圍逐漸縮小, 各工況的差異性逐漸減弱.從圖5 可以觀察到, P =1.1Pc時, 隨溫度的增加, 配位數(shù)曲線的斜率在16.22—5.90區(qū)間變化, 隨著壓力的升高; P =2.0Pc時, 曲線斜率在 14.27—9.25區(qū)間變化.壓力越接近臨界點時, 隨著溫度的變化,密度波動范圍越大, 從而引起配位數(shù)的波動區(qū)間較大, 而在遠離臨界點的壓力下, 密度隨溫度的變化范圍較小, 相應(yīng)的配位數(shù)曲線的斜率范圍較窄.在T/Tpc<1.0的類液區(qū), 壓力增大, 擬臨界點溫度升高, 密度減小, 配位數(shù)斜率呈現(xiàn)減小的變化趨勢;在 T /Tpc≥1.0 的類氣區(qū)表現(xiàn)為相反的變化趨勢.可以得到密度是影響流體配位數(shù)的主要參數(shù), 溫度和壓力也是通過調(diào)整密度的大小間接影響配位數(shù)分布.

圖4 配位數(shù) (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.

根據(jù)配位數(shù)的計算值( Nc)和期望值( Ne)對密度不均勻性進行判斷, 如果直接比較兩個數(shù)字的大小, 這樣的標準過于嚴格, 局部結(jié)構(gòu)的微小波動,即使只有一個分子穿過設(shè)置的半徑邊界, 也會導致密度的分類從平均轉(zhuǎn)變成高或低密度.根據(jù)文獻[30]提到的稍微寬松且能準確判斷的標準, 即選取一個小量的值 δ , 用 Ne± δ 作為參數(shù)度量平均密度的變化, 具體劃分原則如下:

圖5 配位數(shù)曲線斜率Fig.5.The slope of coordination number curve.

計算中允許的波動量為30% Ne, 則 δ 應(yīng)滿足2 δ +1=0.3 Ne, 進一步得到不同參數(shù)下 δ 的具體值,利用劃分原則可以判斷密度分布趨勢, 稱該方法為“30% 方法”.采用該方法, 根據(jù)配位數(shù)得到在P=1.1Pc和 2.0Pc, T =Tpc流體在 xy 平面內(nèi)高/低及平均密度區(qū)分布隨時間的演化如圖6所示.

圖6(a)表示在 0—5000τ的時間范圍內(nèi)P=1.1Pc, T =Tpc時的演化過程, 由圖可知, 該計算工況產(chǎn)生的均值區(qū)面積較小, 低密度區(qū)的位置主要集中在模擬系統(tǒng)的中部, 隨著時間的演化, 密度在不停地波動, 低密度區(qū)呈現(xiàn)分散-聚合-分散的演化規(guī)律, 產(chǎn)生的高密度區(qū)的面積相對較大且位置集中, 演化過程中波動微弱.由于在超臨界區(qū)表面張力的消失, 不存在亞臨界工況下的彎曲界面.形成該現(xiàn)象的原因主要是由于當密度低時, 分子間引力起主導作用, 在溫度的影響下分子處于不斷的相互碰撞中, 能量的增加, 導致高密度區(qū)的形成.此外,該工況擬臨界點位置出現(xiàn)等溫壓縮系數(shù)的極值點,形成較高的密度區(qū)所付出的代價變小, 超臨界流體的關(guān)聯(lián)長度在擬臨界處也存在極大值, 因此會出現(xiàn)大面積的, 相對穩(wěn)定的高密度區(qū).在確定尺寸和分子數(shù)的系統(tǒng)中, 高密度區(qū)的形成必然會引起低密度區(qū)的產(chǎn)生.隨著高密度區(qū)的形成, 分子間的距離縮短, 增大的斥力將部分分子從高密度區(qū)向低密度區(qū)推, 當分子間距離較大時, 引力起主導作用, 會使得分子再次聚集成高密度區(qū), 但此時作用勢效應(yīng)相對等溫壓縮系數(shù)效應(yīng)較小, 因此各區(qū)域所處位置穩(wěn)定, 波動微弱.隨著壓力的升高, 流體的等溫壓縮系數(shù)仍存在極值點, 但是數(shù)值較小, 流體可壓縮性減弱, 較難形成高密度區(qū), 僅在分子間短程勢作用下, 形成由幾個分子組成的且較為分散的高密度區(qū), 形成相對較大的平均值區(qū), 如圖6(b)所示.隨時間的演化, 各區(qū)域不停波動, 存在明顯的嵌套現(xiàn)象, 系統(tǒng)中局部出現(xiàn)類似“花斑”的現(xiàn)象, 這些花斑若隱若現(xiàn)、此起彼伏、互相嵌套的性質(zhì)和相關(guān)教課書[31]中提出的結(jié)論保持一致.

圖6 流體在 xy 平面內(nèi)高/低密度區(qū)分布隨時間的演化 (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 不同壓力下, 擬臨界點溫度下高/低密度區(qū)占比 (a) 高密度區(qū)占比; (b) 低密度區(qū)占比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可以直觀地觀察各工況不同密度區(qū)所占比例隨時間的變化.不同的壓力條件下, 類液和類氣區(qū)的占比一直處于一個波動的狀態(tài), 僅觀察曲線, 發(fā)現(xiàn)波動過程并沒有實際的規(guī)律可循, 整體表現(xiàn)為一個混亂的動態(tài)變化過程.從圖7(a)和圖7(b)分圖中發(fā)現(xiàn)隨著壓力增大類液和類氣區(qū)的占比整體呈現(xiàn)下降的趨勢, 進一步說明隨壓力增大, 形成高密度區(qū)較為困難, 計算壓力距臨界壓力越遠, 流體的均勻性越強.在 P =1.1Pc, T =Tpc的工況時,高密度區(qū)占比約為60%, 低密度區(qū)的占比約為35%,此時均勻區(qū)域的占比最小, 系統(tǒng)的不均勻性最強.

3.2 靜態(tài)結(jié)構(gòu)因子

結(jié)構(gòu)因子主要表征材料對射線的散射能量, 反映材料結(jié)構(gòu)的平均信息, 可以進一步應(yīng)用到流體中, 觀察流體的結(jié)構(gòu)特性.而靜態(tài)結(jié)構(gòu)因子和徑向分布函數(shù)互為傅里葉變換[19,32], 計算式為

其中 n = N/V, V 為系統(tǒng)體積, 由于結(jié)構(gòu)因子是倒易空間, 變量k與距中心分子的距離 rc成反比.

由圖8可知, 各模擬工況下均在較小的k值范圍內(nèi)存在發(fā)散行為, k < 0.5 σ?1, 隨著壓力的增加,這種發(fā)散行為逐漸減弱.在 P =1.1Pc時, 擬臨界點處的發(fā)散行為最為強烈, 流體表現(xiàn)為較強的小角度散射, 曲線在一定范圍內(nèi)存在微小的波動, 而隨著壓力的升高, 這種波動現(xiàn)象消失.在 P =1.3Pc和1.5Pc時, 發(fā)散最強烈的行為出現(xiàn)在偏離擬臨界點的 T =1.01Tpc工況, 與等溫壓縮系數(shù)的最大值點保持一致.在圖8(a)—圖8(c)分圖中得到隨著溫度偏離擬臨界點的距離增大, 發(fā)散性減小, 即各工況均在 T =0.95Tpc和 1.05Tpc時得到靜態(tài)結(jié)構(gòu)因子的最小值.在圖8(d)分中, 由于等溫壓縮系數(shù)在一個較小的水平緩慢增加, 因此靜態(tài)結(jié)構(gòu)因子的發(fā)散特性也呈現(xiàn)微弱增加趨勢.壓力較低時, 靜態(tài)結(jié)構(gòu)因子的發(fā)散在擬臨界點溫度達到最大值, 但隨著壓力的增加, 發(fā)散的極值點工況與等溫壓縮系數(shù)的極值點工況相符合.

3.3 密度時間序列曲線和排列熵

密度隨時間的變化是分子動力學模擬中比較容易得到的數(shù)據(jù), 該曲線給出初步粗略的波動信息.從圖9(a)可以看出, 壓力 P =1.1Pc, 密度時間序列曲線呈現(xiàn)出包含類周期特征的大幅波動, 流體結(jié)構(gòu)具有局部有序特征, 但這種波動存在不規(guī)則運動的相互疊加, 大波套小波的現(xiàn)象, 并不具有嚴格的周期性.隨著壓力的增加, 波動的幅度和類周期性均減弱, 隨機性增強, 此時流體結(jié)構(gòu)具有較強的無序性.

圖8 靜態(tài)結(jié)構(gòu)因子 (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.

為具體描述密度時間序列的復雜程度, 引入排列熵的分析方法, 對一個確定長度的時間序列, 可以根據(jù)自相關(guān)函數(shù)法確定一個延遲時間 τD, 根據(jù)排列熵嵌入維度的要求, 文中選取各工況嵌入維數(shù)m = 5.排列熵作為衡量時間序列曲線復雜度的指標, 越規(guī)則的時間序列, 對應(yīng)的排列熵越小, 越復雜的時間序列, 對應(yīng)的排列熵越大.

根據(jù)文獻[33]提出的計算方法, 對各工況的密度時間序列曲線進行排列熵的計算, 從圖9(b)中可以觀察到在 P =1.1Pc時, 在 T =0.99Tpc和1.0Tpc時排列熵的值較小, 對應(yīng)較規(guī)則的時間序列.隨著壓力的增大, P =1.3Pc和 1.5Pc時, 均在T=1.01Tpc時得到最小的排列熵, 隨著壓力的進一步增加, 排列熵的值呈微弱增大趨勢, 且各工況的值在一個水平線附近波動.產(chǎn)生這種現(xiàn)象的原因主要是分子間作用勢、等溫壓縮系數(shù)和分子熱運動效應(yīng)共同作用的結(jié)果.P =1.1Pc時, 各計算工況的溫度低, 分子熱運動效應(yīng)較弱, 等溫壓縮系數(shù)在擬臨界點處達到極大值, 流體容易被壓縮形成高密度區(qū), 在這兩種效應(yīng)的共同作用下, 排列熵的最小值出現(xiàn)在T=0.99Tpc的工況.隨著壓力的升高( P =1.3Pc和1.5Pc), 擬臨界溫度逐漸增大, 分子的熱運動加劇,導致分子間的作用勢效應(yīng)減小, 此時等溫壓縮系數(shù)效應(yīng)在密度波動中占據(jù)主導地位, 最小排列熵和等溫壓縮系數(shù)極大值點工況一致.當 P =2.0Pc時, 各工況得到的排列熵變化趨勢和等溫壓縮系數(shù)保持一致.

4 結(jié) 論

采用分子動力學方法模擬不同超臨界壓力, 擬臨界溫度附近流體的結(jié)構(gòu)特征, 對流體的徑向分布函數(shù)、配位數(shù)、靜態(tài)結(jié)構(gòu)因子、密度時間序列曲線及排列熵展開研究, 分析壓力和溫度的變化對各參量的影響, 得到如下結(jié)論:

1) 徑向分布函數(shù)在類液區(qū)存在“短程有序、長程無序”的現(xiàn)象, 隨著溫度的升高, 這種有序的現(xiàn)象逐漸消失, 在類氣區(qū)僅在短程區(qū)域存在一個峰值, 而后單調(diào)下降到; 定壓下, 隨著溫度的升高, 第一峰值出現(xiàn)的位置逐漸右移, 但這種變化趨勢隨著壓力偏離臨界點距離的增大而減弱; 此外首次提出, 在 P =1.1Pc和 1.3Pc時, 峰值的最大值均在等溫壓縮極大值點工況出現(xiàn), 但是隨著壓力的增加,該現(xiàn)象不再明顯.

2) 定壓時, 隨著溫度的升高, 配位數(shù)逐漸減小, 壓力的增加導致不同溫度下流體密度的差值減小, 因此引起配位數(shù)的波動范圍縮小, 分布區(qū)域變窄; 采用配位數(shù)標定高/低密度區(qū)的標準, 在P=1.1Pc, T =Tpc時可以得到大面積、分布集中且波動較小的高/低密度區(qū), 此時均值區(qū)域占比較小; 隨著壓力升高, 均值區(qū)域逐漸增加, 高/低密度區(qū)逐漸減小, 僅有幾個分子大小, 且相互嵌套, 并隨時間有較明顯的波動, 說明隨著壓力的增大, 流體的均勻性逐漸增強.

3) 靜態(tài)結(jié)構(gòu)因子是通過散射效應(yīng)觀察流體的內(nèi)部結(jié)構(gòu), 各模擬工況均在較小的k值內(nèi)存在發(fā)散行為, 隨著壓力的增加, 靜態(tài)結(jié)構(gòu)因子曲線的發(fā)散值呈迅速下降的趨勢, 研究發(fā)現(xiàn)定壓時靜態(tài)結(jié)構(gòu)因子最大值和等溫壓縮系數(shù)極值點工況保持一致.

4) 密度時間序列曲線可以初步得到中心切片密度隨壓力的增加, 波動幅度逐漸減弱的規(guī)律, 和前文壓力增大流體均勻性增強的結(jié)論符合.排列熵的變化規(guī)律可以分為三類: 一是在 P =1.1Pc的低壓下, 排列熵在小于擬臨界溫度的位置得到最大值;二是在 P =1.3Pc和 1.5Pc的中壓下, 排列熵的最小值點和等溫壓縮系數(shù)極值點一致; 三是在P =2.0Pc的高壓下, 排列熵的值在一個水平線附近波動, 整體變化趨勢較為平緩, 流體均勻性增強.探討微尺度下超臨界流體的結(jié)構(gòu)特征有助于全面了解超臨界流體特性, 為超臨界流體的應(yīng)用提供有力支撐.