超臨界水中氫鍵特性及其空間分布的密度相關性

曾冬雷, 馮 飆, 宋嘉文, 范利武,2

(1. 浙江大學 熱工與動力系統(tǒng)研究所, 浙江 杭州 310027;2. 浙江大學 能源清潔利用國家重點實驗室, 浙江 杭州 310027)

1 前 言

超臨界流體是溫度高于臨界溫度Tc與壓力高于其臨界壓力pc的特殊流體，由于其普遍具有較高的氧化性、溶解度和擴散系數(shù)，因此常被用做氧化劑[1-2]或催化劑[3]等反應介質。近年來，隨著對超臨界流體物性研究的不斷深入，其密度的局部不均勻性受到了研究人員的廣泛關注，并采用分子動力學方法進行了研究。分子動力學方法可從微觀機理角度解釋宏觀現(xiàn)象并且不受實驗技術限制，可模擬極端條件下的體系并計算得到相關物理量[4]，是研究以超臨界水(臨界溫度與壓力分別為647.15 K 與22.1 MPa)為代表的超臨界流體物性的有力工具。

SKARMOUTSOS 等[5]發(fā)現(xiàn)溫度為 666 K 的超臨界水(模擬的水分子數(shù)量為 500 個，密度范圍為 0.2 ρc～2 ρc，ρc為水的臨界密度)存在局部密度與結構的不均勻性，而且含氫鍵流體的局部密度不均勻性更大[5-6]。氫鍵的鍵能介于共價鍵和范德華力之間，考慮分子間氫鍵的平均強度，形成一個氫鍵則增強了其周圍氫鍵的協(xié)同作用，且氫鍵的協(xié)同作用可使水分子以大的團簇結構存在[7]。MA 等[8]從氫鍵數(shù)量與體系能量波動來解釋這種不均勻性，并發(fā)現(xiàn)溫度在666 K 的條件下，密度在低于0.19 ρc時能量波動急劇增加，此時四面體序參數(shù)有最小值，體系短程有序被破壞。SKARMOUTSOS 等[9]發(fā)現(xiàn)當溫度為666 K 時，周圍有1 個和2 個氫鍵的水分子對超臨界水的局部密度不均勻性的貢獻最大，并從三角形和四面體序參數(shù)角度揭示了局部有序性與氫鍵的關系。SAHLE 等[10]則通過從頭分子動力學模擬和密度泛函理論方法，發(fā)現(xiàn)在溫度低于873 K 與壓力低于132 MPa 的超臨界水中仍然存在著氫鍵網絡。

由于密度不均勻性可能引起超臨界水中氫鍵特性的變化，并且氫鍵作為貢獻較大的分子間相互作用，其數(shù)量、鍵能大小及空間分布會對超臨界水的物性造成重要影響。因此對氫鍵特性及其空間分布隨密度變化的研究尤為重要。但是目前文獻中關于超臨界水密度對氫鍵鍵能及其空間分布影響的研究較為缺乏。因此，本文參照現(xiàn)有文獻的研究結果，對溫度為近臨界溫度666 K 下不同密度的超臨界水體系進行分子動力學模擬，并增加更高溫度的700 與800 K 這2 個工況作為對照組，研究不同溫度下超臨界水的密度對其氫鍵數(shù)量、鍵能以及空間分布規(guī)律的影響，將有助于進一步理解超臨界水密度不均勻性的微觀機制。

2 分子動力學模擬

2.1 模擬細節(jié)

本研究采用經典的TIP4P/2005[11]作為水分子模型。如圖1 所示為超臨界水的初始分子建模示意圖，其中紅球表示氧原子，綠球表示氫原子，立方體的邊長為L。所模擬的超臨界水體系含1 000 個水分子，不同密度下模型的L 取值如表1 所示(ρ/ρc為超臨界水密度與水臨界密度的比值)。

分子動力學模擬使用的平臺是LAMMPS，采用的力場為TIP4P 力場，并設置三維周期性邊界條件(PPP)。溫度與壓力的控制采用Nose-Hoover 耦合方法，并采用PPPM 算法計算靜電相互作用力。初始分子速度由高斯隨機分布決定。庫倫截斷半徑與范德華力截斷半徑均為1 nm。根據(jù)相關文獻中對水分子進行模擬的經驗，模擬的時間步長設為10-15s，每隔10-12s 保存坐標數(shù)據(jù)。在正則系綜(即NVT 系綜)下進行模擬，系統(tǒng)先運行2×10-9s 達到平衡態(tài)，再計算10-9s 進行數(shù)據(jù)分析。如上所述，為揭示不同溫度的影響，系統(tǒng)的溫度設定為666、700、800 K。

圖1 超臨界水的初始分子建模Fig.1 Initial molecular model of supercritical water

表1 超臨界水模型的邊長取值Table 1 Side length values of the supercritical water model

2.2 模型檢驗

為了對模擬體系的分子數(shù)量無關性進行檢驗，首先計算了溫度為666 K 時分子數(shù)分別為500、1 000與5 000 的3 個超臨界水體系中的平均氫鍵數(shù)隨密度的變化關系。如圖2 所示，3 個體系模擬得到的平均氫鍵數(shù)(＜ nHB＞permolecule)隨密度比值的變化規(guī)律完全一致，體系分子數(shù)的變化對模擬結果影響不大，相對偏差在20%以內。因此，為了相對節(jié)約計算資源，本研究的后續(xù)模擬中均選取水分子數(shù)為1 000 來構建模擬體系。

對模擬得到的超臨界水的自擴散系數(shù)與文獻中的實驗值[12]以及前人的分子動力學模擬結果(采用 256個水分子，溫度為666 K)[8]進行對照，以檢驗本文所建立的超臨界水模型的可靠性。采用Einstein 法[13]計算超臨界水的自擴散系數(shù)。該方法的基本原理是通過統(tǒng)計粒子均方根位移與運動時間的關系，當時間足夠長二者呈線性關系時取其斜率的1/6 即為自擴散系數(shù)，其公式如下所示：

圖2 不同水分子數(shù)超臨界水體系中平均氫鍵數(shù)隨密度的變化Fig.2 Variation of average number of hydrogen bonds as a function of density in supercritical water systems with different water molecule numbers

式中：N 為粒子數(shù)目，t 為運動時間，ri(t)為粒子瞬時位置，ri(0)為粒子初始位置。

實驗值由LAMB 等[12]的實驗數(shù)據(jù)插值計算得到。如圖3 所示，本文的模擬值與實驗值和前人的模擬值均吻合較好，總體偏差率低于15%。在保持溫度不變的條件下，超臨界水的自擴散系數(shù)D 隨密度比值升高而減小至趨于平緩。這是由于密度增大，分子間的碰撞概率增大，且氫鍵數(shù)量增加，進一步阻礙了分子擴散，傳質能力減弱，從而表現(xiàn)為自擴散系數(shù)減小。在ρ ＞ 1.2 ρc時自擴散系數(shù)的變化低于 10%，與 MA 等[8]在ρ ＞ 1.25 ρc時發(fā)現(xiàn)的密度對自擴散系數(shù)影響很小的結果基本一致。800 K 超臨界水自擴散系數(shù)最高，666 K 最低，且隨著溫度的降低自擴散系數(shù)也逐漸減小，是由于溫度升高促進了水分子的熱運動，且氫鍵作用減弱，水分子運動的空間位阻效應減弱[14]，傳質能力增強。

2.3 徑向分布函數(shù)

如圖4 所示為徑向分布示意圖，通常用來描述局部粒子的分布，其表達式如下所示：

圖3 超臨界水自擴散系數(shù)隨密度的變化Fig.3 Variation of self-diffusion coefficient of supercritical water as a function of density

圖4 徑向分布示意圖Fig.4 Schematic diagram of radial distribution

式中：dnr為半徑從r 到r+dr 球殼內圍繞每個A 粒子周圍B 粒子的數(shù)目，ρB為B 粒子的數(shù)密度，r 為距中心原子的半徑，dr 為球殼厚度。徑向分布函數(shù)實際上是A 粒子周圍B 粒子出現(xiàn)的概率密度函數(shù)，可通過其觀測粒子的空間分布。由于每個水分子中只包含 1 個氧原子，因此可以通過觀察氧氧徑向分布函數(shù)來描述水分子的分布狀況，即從 gOO(r)峰值出現(xiàn)的位置，高度及曲線的走勢來分析水分子的分布狀況。因為水分子的分布狀況的變化會影響氫鍵的分布狀況，所以分析水分子分布對氫鍵分布的研究具有重要意義。

2.4 氫鍵特性表征

本文使用幾何標準定義氫鍵[6]，從而對計算輸出的坐標文件進行編程統(tǒng)計氫鍵的數(shù)量。如圖5 所示，供體與受體氧原子距離rOO≤0.36 nm，供體與受體氧原子與形成氫鍵的氫原子夾角∠HOO≤30°且受體氧原子與氫原子距離rOH≤0.24 nm 可視為形成氫鍵。

圖5 氫鍵幾何模型Fig.5 Geometric model of hydrogen bonds

氫鍵鍵能僅由3 個原子(OD-H...OA)的短程庫侖作用決定：

3 模擬結果及討論

3.1 水分子的空間分布

由于每個水分子只含有一個氧原子，因此可以通過計算氧氧徑向分布來分析水分子的空間分布情況。不同密度超臨界水中的氧氧徑向分布模擬結果如圖6 所示。

圖6 超臨界水中的氧氧徑向分布Fig.6 Radial distribution of oxygen-oxygen in supercritical water

根據(jù)NVT 系綜下分子運動軌跡文件，計算了在666，700與800 K 溫度下的不同密度的氧氧原子的徑向分布函數(shù)。由圖6可知，超臨界水gOO(r)的第1 個峰均位于0.29 nm，說明超臨界條件下水分子之間出現(xiàn)概率最高為0.29 nm。第2 峰基本消失，說明在超臨界條件下氫鍵作用變得很弱。在恒溫條件下，gOO(r)的第1 峰高度具有密度依賴性，且隨著密度的減小，第1 峰高度增大，說明超臨界水短程有序性隨密度減小而增大，在較低密度的條件下會形成局部團聚現(xiàn)象。水分子的局部團聚會導致氫鍵在空間分布上呈現(xiàn)局部團聚的現(xiàn)象。但隨著溫度升高，由于分子熱運動增強，這種密度依賴性減弱，且第1 峰高度減小，這與SOPER 等[15]的研究結論一致。

3.2 氫鍵數(shù)量與平均鍵能

水分子平均氫鍵數(shù)(＜ nHB＞permolecule)與密度比值的關系如圖7 所示。值得注意的是，圖中給出了在若干次模擬中所得到數(shù)據(jù)點的誤差波動情況以說明模擬數(shù)據(jù)的可重復性。由于使用了NVT 系綜，因此在模擬過程中存在體系溫度的隨機波動現(xiàn)象，從而導致模擬結果也會出現(xiàn)一定的隨機波動。在溫度不變的情況下，水分子平均氫鍵數(shù)隨密度增加而增加，這是由于水分子密度增大，結構更緊湊，水分子之間更容易形成氫鍵。但當 ρ ＞1.4ρc時平均氫鍵數(shù)趨于平緩，密度的影響減小。超臨界水在666 K 下平均氫鍵數(shù)最少，800 K時最多。這是由于水分子的熱運動增強導致更難達到成鍵的幾何條件。因此，隨著溫度的升高，氫鍵的作用減弱。在666 K 溫度下，平均氫鍵數(shù)在密度為2.0ρc時的平均氫鍵數(shù)較高，可能是因為在較低溫度下氫鍵的密度依賴性更大，高密度時更容易出現(xiàn)團聚現(xiàn)象。

圖7 超臨界水中氫鍵數(shù)量隨密度的變化Fig.7 Variation of hydrogen bond number as a function of density in supercritical water

圖8 超臨界水中平均氫鍵鍵能隨密度的變化Fig.8 Variation of average hydrogen bond energy as a function of density in supercritical water

平均氫鍵鍵能EHB與密度比值的關系如圖8 所示。在溫度一定的條件下，超臨界水平均氫鍵鍵能隨著密度比值增大有上升的趨勢，其中溫度為666 K 時平均氫鍵鍵能在2.0 ρc相比0.4 ρc增加了7.0%，溫度為700 K 時增加了9.1%，溫度為800 K 時增加了12.3%。氫鍵鍵能增長幅度隨溫度升高而增加，總體增長幅度較為平緩。這是因為密度增大導致水分子之間距離更近，形成的氫鍵更穩(wěn)定，鍵能更大。在溫度為800 K 的條件下，超臨界水平均氫鍵鍵能則波動最大，可能是由于水分子熱運動較強，使得原子之間的距離有很大的隨機性；在溫度為666 K 的條件下，平均氫鍵鍵能總體最高，且隨和溫度的上升而下降，這是因為溫度更低的條件下分子間距離更小，導致氫鍵鍵能更大。

3.3 氫鍵的空間分布

圖9 不同密度(0.4 ρc，1.2 ρc，2.0 ρc)超臨界水中氫鍵的空間分布Fig.9 Spatial distribution of hydrogen bonds in supercritical water having different densities (0.4 ρc，1.2 ρc，2.0 ρc)

如圖9 所示以受體氧原子坐標為氫鍵坐標，氫鍵鍵能大小表現(xiàn)在圓球的大小上。溫度不變時，隨著密度的增大，氫鍵成鍵在空間分布上更均勻，團聚現(xiàn)象(圖中黑色虛線所示)越少，密度增大與氫鍵更均勻的分布使得水分子的傳質能力降低，表現(xiàn)在自擴散系數(shù)的減小。高密度超臨界水體系中心處平均氫鍵鍵能大小更高(圓球更大)，成鍵更穩(wěn)定。在密度不變的條件下，氫鍵空間分布隨著溫度的升高而變得更加均勻，即更難出現(xiàn)團聚現(xiàn)象。同時，在更高溫度的情況下，氫鍵數(shù)目更低，氫鍵鍵能大小差異更小。

4 結 論

本文采用分子動力學模擬方法分析了溫度為666，700 與800 K 的超臨界水在不同密度下(0.4 ρc～2 ρc)的氫鍵特性及其空間分布與密度的關系，得出以下主要結論：

(1) 超臨界水中的氫鍵數(shù)量隨其密度的增大而增加，但在ρ ＞ 1.4ρc時氫鍵數(shù)量的增長較為平緩；在較高的溫度下氫鍵數(shù)量也較少。

(2) 超臨界水中氫鍵的平均鍵能隨其密度的增大而增大。當溫度升高時，平均氫鍵鍵能變小，但其隨密度變化的波動性變得更大。

(3) 超臨界水中氫鍵的成鍵空間均勻性隨密度的增大而增加，且隨溫度的升高而增加。密度較低時，超臨界水更易形成局部團聚現(xiàn)象，導致氫鍵在空間分布上也呈現(xiàn)局部團聚；而密度較高時，超臨界水體系在幾何中心處的氫鍵成鍵更穩(wěn)定。