王艷 徐進(jìn)良2)? 李文 劉歡

1) (華北電力大學(xué), 低品位能源多相流與傳熱北京市重點(diǎn)實(shí)驗(yàn)室, 北京 102206)

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

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

1 引 言

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

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

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

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

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

超臨界流體氬的臨界點(diǎn)溫度(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個工況進(jìn)行模擬分析.

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

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

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

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

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

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

3 結(jié)果與討論

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

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

其中 ρ 為系統(tǒng)的密度, N為分子的總數(shù)目, D為計(jì)算的總時間 (步數(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得知, 在各計(jì)算工況下, 溫度低于擬臨界點(diǎn)溫度時, 徑向分布函數(shù)會出現(xiàn)高低不同的峰谷, 在1附近產(chǎn)生振蕩, 表現(xià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é)果.氬的局部聚集必然會導(dǎo)致體系內(nèi)部出現(xiàn)“空隙”, 壓力的增大將這些“空隙”壓縮, 減小了局部聚集氬相對密度, 因而會出現(xiàn)峰值和峰寬度都減小的現(xiàn)象, 和現(xiàn)有文獻(xiàn)[27]分析保持一致.同時也進(jìn)一步證明超臨界流體的高壓縮性.在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)在最低溫度工況下, 這主要由于在擬臨界點(diǎn)附近超臨界流體存在很大的壓縮性,促使流體聚集產(chǎn)生高密度區(qū), 但是隨著壓力的增加, 流體的壓縮性迅速降低, 溫度升高帶來的分子熱運(yùn)動逐漸加強(qiáng), 使得流體的聚集能力逐漸減弱.

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

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

有研究表明配位數(shù)不僅與結(jié)構(gòu)、密度有關(guān), 而且是溫度的函數(shù)[29].在相同的壓力下, 隨著溫度的升高, 流體的密度減小, 此時分子間的距離增大,分子間引力對配位數(shù)的影響越來越大, 但是溫度的升高, 促使分子的熱運(yùn)動加劇, 在引力和熱運(yùn)動的共同作用下, 發(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ū)間變化.壓力越接近臨界點(diǎn)時, 隨著溫度的變化,密度波動范圍越大, 從而引起配位數(shù)的波動區(qū)間較大, 而在遠(yuǎn)離臨界點(diǎn)的壓力下, 密度隨溫度的變化范圍較小, 相應(yīng)的配位數(shù)曲線的斜率范圍較窄.在T/Tpc<1.0的類液區(qū), 壓力增大, 擬臨界點(diǎn)溫度升高, 密度減小, 配位數(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ù)的計(jì)算值( Nc)和期望值( Ne)對密度不均勻性進(jìn)行判斷, 如果直接比較兩個數(shù)字的大小, 這樣的標(biāo)準(zhǔn)過于嚴(yán)格, 局部結(jié)構(gòu)的微小波動,即使只有一個分子穿過設(shè)置的半徑邊界, 也會導(dǎo)致密度的分類從平均轉(zhuǎn)變成高或低密度.根據(jù)文獻(xiàn)[30]提到的稍微寬松且能準(zhǔn)確判斷的標(biāo)準(zhǔn), 即選取一個小量的值 δ , 用 Ne± δ 作為參數(shù)度量平均密度的變化, 具體劃分原則如下:

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

計(jì)算中允許的波動量為30% Ne, 則 δ 應(yīng)滿足2 δ +1=0.3 Ne, 進(jìn)一步得到不同參數(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時的演化過程, 由圖可知, 該計(jì)算工況產(chǎn)生的均值區(qū)面積較小, 低密度區(qū)的位置主要集中在模擬系統(tǒng)的中部, 隨著時間的演化, 密度在不停地波動, 低密度區(qū)呈現(xiàn)分散-聚合-分散的演化規(guī)律, 產(chǎn)生的高密度區(qū)的面積相對較大且位置集中, 演化過程中波動微弱.由于在超臨界區(qū)表面張力的消失, 不存在亞臨界工況下的彎曲界面.形成該現(xiàn)象的原因主要是由于當(dāng)密度低時, 分子間引力起主導(dǎo)作用, 在溫度的影響下分子處于不斷的相互碰撞中, 能量的增加, 導(dǎo)致高密度區(qū)的形成.此外,該工況擬臨界點(diǎn)位置出現(xiàn)等溫壓縮系數(shù)的極值點(diǎn),形成較高的密度區(qū)所付出的代價變小, 超臨界流體的關(guān)聯(lián)長度在擬臨界處也存在極大值, 因此會出現(xiàn)大面積的, 相對穩(wěn)定的高密度區(qū).在確定尺寸和分子數(shù)的系統(tǒng)中, 高密度區(qū)的形成必然會引起低密度區(qū)的產(chǎn)生.隨著高密度區(qū)的形成, 分子間的距離縮短, 增大的斥力將部分分子從高密度區(qū)向低密度區(qū)推, 當(dāng)分子間距離較大時, 引力起主導(dǎo)作用, 會使得分子再次聚集成高密度區(qū), 但此時作用勢效應(yīng)相對等溫壓縮系數(shù)效應(yīng)較小, 因此各區(qū)域所處位置穩(wěn)定, 波動微弱.隨著壓力的升高, 流體的等溫壓縮系數(shù)仍存在極值點(diǎn), 但是數(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 不同壓力下, 擬臨界點(diǎn)溫度下高/低密度區(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)波動過程并沒有實(shí)際的規(guī)律可循, 整體表現(xiàn)為一個混亂的動態(tài)變化過程.從圖7(a)和圖7(b)分圖中發(fā)現(xiàn)隨著壓力增大類液和類氣區(qū)的占比整體呈現(xiàn)下降的趨勢, 進(jìn)一步說明隨壓力增大, 形成高密度區(qū)較為困難, 計(jì)算壓力距臨界壓力越遠(yuǎn), 流體的均勻性越強(qiáng).在 P =1.1Pc, T =Tpc的工況時,高密度區(qū)占比約為60%, 低密度區(qū)的占比約為35%,此時均勻區(qū)域的占比最小, 系統(tǒng)的不均勻性最強(qiáng).

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

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

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

由圖8可知, 各模擬工況下均在較小的k值范圍內(nèi)存在發(fā)散行為, k < 0.5 σ?1, 隨著壓力的增加,這種發(fā)散行為逐漸減弱.在 P =1.1Pc時, 擬臨界點(diǎn)處的發(fā)散行為最為強(qiáng)烈, 流體表現(xiàn)為較強(qiáng)的小角度散射, 曲線在一定范圍內(nèi)存在微小的波動, 而隨著壓力的升高, 這種波動現(xiàn)象消失.在 P =1.3Pc和1.5Pc時, 發(fā)散最強(qiáng)烈的行為出現(xiàn)在偏離擬臨界點(diǎn)的 T =1.01Tpc工況, 與等溫壓縮系數(shù)的最大值點(diǎn)保持一致.在圖8(a)—圖8(c)分圖中得到隨著溫度偏離擬臨界點(diǎn)的距離增大, 發(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ā)散在擬臨界點(diǎn)溫度達(dá)到最大值, 但隨著壓力的增加, 發(fā)散的極值點(diǎn)工況與等溫壓縮系數(shù)的極值點(diǎn)工況相符合.

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

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

圖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.

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

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

4 結(jié) 論

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

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

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

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ù)極值點(diǎn)工況保持一致.

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