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

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

1 前 言

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

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

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

2 分子動力學(xué)模擬

2.1 模擬細(xì)節(jié)

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

分子動力學(xué)模擬使用的平臺是LAMMPS，采用的力場為TIP4P 力場，并設(shè)置三維周期性邊界條件(PPP)。溫度與壓力的控制采用Nose-Hoover 耦合方法，并采用PPPM 算法計算靜電相互作用力。初始分子速度由高斯隨機(jī)分布決定。庫倫截斷半徑與范德華力截斷半徑均為1 nm。根據(jù)相關(guān)文獻(xiàn)中對水分子進(jìn)行模擬的經(jīng)驗，模擬的時間步長設(shè)為10-15s，每隔10-12s 保存坐標(biāo)數(shù)據(jù)。在正則系綜(即NVT 系綜)下進(jìn)行模擬，系統(tǒng)先運(yùn)行2×10-9s 達(dá)到平衡態(tài)，再計算10-9s 進(jìn)行數(shù)據(jù)分析。如上所述，為揭示不同溫度的影響，系統(tǒng)的溫度設(shè)定為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ù)量無關(guān)性進(jìn)行檢驗，首先計算了溫度為666 K 時分子數(shù)分別為500、1 000與5 000 的3 個超臨界水體系中的平均氫鍵數(shù)隨密度的變化關(guān)系。如圖2 所示，3 個體系模擬得到的平均氫鍵數(shù)(＜ nHB＞permolecule)隨密度比值的變化規(guī)律完全一致，體系分子數(shù)的變化對模擬結(jié)果影響不大，相對偏差在20%以內(nèi)。因此，為了相對節(jié)約計算資源，本研究的后續(xù)模擬中均選取水分子數(shù)為1 000 來構(gòu)建模擬體系。

對模擬得到的超臨界水的自擴(kuò)散系數(shù)與文獻(xiàn)中的實驗值[12]以及前人的分子動力學(xué)模擬結(jié)果(采用 256個水分子，溫度為666 K)[8]進(jìn)行對照，以檢驗本文所建立的超臨界水模型的可靠性。采用Einstein 法[13]計算超臨界水的自擴(kuò)散系數(shù)。該方法的基本原理是通過統(tǒng)計粒子均方根位移與運(yùn)動時間的關(guān)系，當(dāng)時間足夠長二者呈線性關(guān)系時取其斜率的1/6 即為自擴(kuò)散系數(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 為運(yùn)動時間，ri(t)為粒子瞬時位置，ri(0)為粒子初始位置。

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

2.3 徑向分布函數(shù)

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

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

2.4 氫鍵特性表征

本文使用幾何標(biāo)準(zhǔn)定義氫鍵[6]，從而對計算輸出的坐標(biāo)文件進(jìn)行編程統(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 模擬結(jié)果及討論

3.1 水分子的空間分布

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

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

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

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

水分子平均氫鍵數(shù)(＜ nHB＞permolecule)與密度比值的關(guān)系如圖7 所示。值得注意的是，圖中給出了在若干次模擬中所得到數(shù)據(jù)點的誤差波動情況以說明模擬數(shù)據(jù)的可重復(fù)性。由于使用了NVT 系綜，因此在模擬過程中存在體系溫度的隨機(jī)波動現(xiàn)象，從而導(dǎo)致模擬結(jié)果也會出現(xiàn)一定的隨機(jī)波動。在溫度不變的情況下，水分子平均氫鍵數(shù)隨密度增加而增加，這是由于水分子密度增大，結(jié)構(gòu)更緊湊，水分子之間更容易形成氫鍵。但當(dāng) ρ ＞1.4ρc時平均氫鍵數(shù)趨于平緩，密度的影響減小。超臨界水在666 K 下平均氫鍵數(shù)最少，800 K時最多。這是由于水分子的熱運(yùn)動增強(qiáng)導(dǎo)致更難達(dá)到成鍵的幾何條件。因此，隨著溫度的升高，氫鍵的作用減弱。在666 K 溫度下，平均氫鍵數(shù)在密度為2.0ρc時的平均氫鍵數(shù)較高，可能是因為在較低溫度下氫鍵的密度依賴性更大，高密度時更容易出現(xiàn)團(tuá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與密度比值的關(guān)系如圖8 所示。在溫度一定的條件下，超臨界水平均氫鍵鍵能隨著密度比值增大有上升的趨勢，其中溫度為666 K 時平均氫鍵鍵能在2.0 ρc相比0.4 ρc增加了7.0%，溫度為700 K 時增加了9.1%，溫度為800 K 時增加了12.3%。氫鍵鍵能增長幅度隨溫度升高而增加，總體增長幅度較為平緩。這是因為密度增大導(dǎo)致水分子之間距離更近，形成的氫鍵更穩(wěn)定，鍵能更大。在溫度為800 K 的條件下，超臨界水平均氫鍵鍵能則波動最大，可能是由于水分子熱運(yùn)動較強(qiáng)，使得原子之間的距離有很大的隨機(jī)性；在溫度為666 K 的條件下，平均氫鍵鍵能總體最高，且隨和溫度的上升而下降，這是因為溫度更低的條件下分子間距離更小，導(dǎo)致氫鍵鍵能更大。

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

4 結(jié) 論

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

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

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

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