超臨界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的高壓下, 排列熵的值在一個水平線附近波動, 整體變化趨勢較為平緩, 流體均勻性增強.探討微尺度下超臨界流體的結構特征有助于全面了解超臨界流體特性, 為超臨界流體的應用提供有力支撐.