超臨界態(tài)二氧化碳流體熱力學(xué)性質(zhì)的分子動(dòng)力學(xué)模擬

劉麗麗 代威 楊智 廖傳軍

(1 西安航天神舟建筑設(shè)計(jì)院有限公司北京分公司 北京 100040)

(2 北京宇航系統(tǒng)工程研究所 北京 100076)

(3 廣東工業(yè)大學(xué)材料與能源學(xué)院 廣州 510006)

1 引言

二氧化碳化學(xué)性質(zhì)穩(wěn)定,無毒,不易燃易爆,臭氧消耗和氣候變暖潛值較低,天然綠色環(huán)保,臨界狀態(tài)易于實(shí)現(xiàn)。在臨界點(diǎn)以上,二氧化碳具有類似于液體的密度,氣體的粘度和擴(kuò)散系數(shù),溶解能力強(qiáng),表面張力為零,具有良好的傳熱傳質(zhì)特性[1]。目前,CO2在超臨界萃取、制冷熱泵及熱力發(fā)電系統(tǒng)、頁(yè)巖壓裂及強(qiáng)化石油開采和環(huán)境保護(hù)(如CO2捕集和封存)等諸多領(lǐng)域已有廣泛應(yīng)用。然而,超臨界態(tài)流體熱物性的準(zhǔn)確獲取是其得以廣泛應(yīng)用的重要前提。

傳統(tǒng)上,流體熱物性可通過實(shí)驗(yàn)測(cè)量以及狀態(tài)方程的關(guān)聯(lián)計(jì)算來獲取。然而,為獲取超臨界態(tài)熱物性數(shù)據(jù),實(shí)驗(yàn)測(cè)量方法則需針對(duì)測(cè)量裝置的耐高溫高壓性能進(jìn)行特殊設(shè)計(jì),同時(shí)還需保證較高的測(cè)量精度并測(cè)量足夠多的數(shù)據(jù)點(diǎn),其代價(jià)通常極為昂貴且技術(shù)上難以實(shí)現(xiàn)。對(duì)狀態(tài)方程而言,其通常對(duì)實(shí)驗(yàn)數(shù)據(jù)依賴性較強(qiáng),如專用型狀態(tài)方程[2]的確定則需較多的參數(shù)自由度及大量的實(shí)驗(yàn)數(shù)據(jù)樣本,其適用范圍也嚴(yán)格受限于實(shí)驗(yàn)數(shù)據(jù)樣本所處的熱力狀態(tài)區(qū)間,外推預(yù)測(cè)能力較弱,而簡(jiǎn)單的立方型狀態(tài)方程通常外推預(yù)測(cè)精度相對(duì)不足。

本質(zhì)上,流體的宏觀熱物性由流體分子的微觀相互作用(分子力場(chǎng))所決定。作為熱物性獲取的一種有效方法,基于力場(chǎng)的分子模擬具有堅(jiān)實(shí)的統(tǒng)計(jì)熱力學(xué)理論基礎(chǔ)。近年來,隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,分子模擬在流體熱物性研究領(lǐng)域扮演越來越重要的角色,其不僅可提供分子體系微觀結(jié)構(gòu)的可視化研究,實(shí)現(xiàn)流體宏觀熱物性的微觀機(jī)理解釋,也可作為實(shí)驗(yàn)測(cè)量方法的一種有效補(bǔ)充,獲得具有定量精度的內(nèi)插及外推物性數(shù)據(jù)、獲得實(shí)驗(yàn)無法測(cè)量或難以發(fā)現(xiàn)的結(jié)果(如極端的高溫高壓條件)、預(yù)測(cè)現(xiàn)有或新材料的性質(zhì)及發(fā)展新理論。分子模擬在流體熱物性研究領(lǐng)域已被公認(rèn)為聯(lián)系實(shí)驗(yàn)數(shù)據(jù)和理論模型的重要橋梁。

分子力場(chǎng)是分子模擬的關(guān)鍵,目前國(guó)內(nèi)外學(xué)者針對(duì)CO2流體已開發(fā)出多種不同的分子力場(chǎng),如簡(jiǎn)單的粗?；瘑瘟Ｗ幽Ｐ蚚3]以及高分辨率的全原子模型[4-6]均得到了廣泛應(yīng)用。盡管CO2分子力場(chǎng)存在諸多不同的描述形式,但大多數(shù)力場(chǎng)中描述分子間范德華相互作用的參數(shù)均通過擬合相平衡實(shí)驗(yàn)數(shù)據(jù)來獲取,以實(shí)現(xiàn)其對(duì)亞臨界態(tài)熱物性較好的預(yù)測(cè)精度。然而據(jù)文獻(xiàn)調(diào)研發(fā)現(xiàn),基于氣液相平衡數(shù)據(jù)擬合的半經(jīng)驗(yàn)力場(chǎng)模型對(duì)超臨界態(tài)熱物性的外延預(yù)測(cè)性能仍少有研究。因此,采用分子動(dòng)力學(xué)模擬方法,針對(duì)幾種典型的粗?；腿恿?chǎng)模型,如SAFT-γ[3]、Zhang[4]、EMP2[5]、TraPPE[6]等,開展超臨界態(tài)CO2流體熱力學(xué)特性的預(yù)測(cè)研究,以充分評(píng)估不同力場(chǎng)對(duì)CO2流體超臨界態(tài)熱力學(xué)性質(zhì)的外延預(yù)測(cè)性能,從而為極端溫度和壓力工況下CO2流體熱物性的獲取提供一種有效的替代方案。

2 模擬方法

2.1 力場(chǎng)模型

本研究中,不同力場(chǎng)的范德華相互作用均可采用一般化Mie 勢(shì)能函數(shù)[7]進(jìn)行描述,其函數(shù)形式為:

對(duì)于帶電粒子模型,分子間還需考慮長(zhǎng)程靜電相互作用,如式(3),其中qi和qj為位于原子中心的點(diǎn)電荷,ε0為真空介電常數(shù)。對(duì)于全原子力場(chǎng)模型,體系勢(shì)能還包含分子內(nèi)原子間相互作用,也即鍵伸縮項(xiàng)和鍵角彎曲項(xiàng),如式(4)。其中,kr、kθ分別為鍵伸縮和鍵角彎曲力常數(shù),r0、θ0分別為平衡鍵長(zhǎng)和鍵角。主要考慮了4 種CO2分子力場(chǎng),其中TraPPE、EPM2和Zhang 模型為三點(diǎn)硬球模型,SAFT-γ 模型為單點(diǎn)粗?；Ｐ汀Ｏ啾扔谌c(diǎn)模型,單點(diǎn)模型通常具有更簡(jiǎn)單的結(jié)構(gòu)和更高的計(jì)算效率,在相同的計(jì)算條件下,可用于更大尺度流體行為和現(xiàn)象的模擬。4 種CO2力場(chǎng)參數(shù)如表1 所示;TraPPE 和Zhang 模型采用Lorentz-Berthelot 混合規(guī)則描述不同類原子間相互作用,而EPM2 力場(chǎng)采用幾何平均數(shù)規(guī)則。

表1 不同CO2 力場(chǎng)模型對(duì)應(yīng)的相互作用參數(shù)Table 1 Interaction parameters of different CO2 field models

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

所有模擬均在開源的LAMMPS 平臺(tái)下進(jìn)行,對(duì)于單點(diǎn)粗粒化力場(chǎng)模型,模擬體系包含1 000 個(gè)粒子,對(duì)于全原子力場(chǎng)模型,模擬體系由512 分子填充。模擬均在恒溫恒壓(NPT)系綜下進(jìn)行,且采用周期性邊界條件。模擬的時(shí)間步長(zhǎng)設(shè)置為1 飛秒,同時(shí)采用Nosé-Hoover 方法對(duì)模擬體系進(jìn)行溫度和壓力控制。在數(shù)據(jù)采集階段,采用PPPM 方法計(jì)算長(zhǎng)程庫(kù)倫相互作用,且將靜電力計(jì)算的相對(duì)偏差控制在1 ×10-4以內(nèi)。每個(gè)工況點(diǎn)下,模擬的總步數(shù)設(shè)置為750 萬步,其中前500 萬步用于系統(tǒng)的弛豫平衡,后250 萬步用于數(shù)據(jù)采用及熱力學(xué)性質(zhì)的統(tǒng)計(jì)平均。

2.3 熱物性參數(shù)計(jì)算

在二氧化碳的熱物性參數(shù)計(jì)算中,流體的密度ρ可直接通過系統(tǒng)體積的統(tǒng)計(jì)平均值來計(jì)算,也即:

式中:N為模擬體系的分子數(shù),M為分子摩爾質(zhì)量,V為模擬系統(tǒng)體積,Na為阿伏伽德羅常數(shù)。

通常,流體的比定壓熱容Cp由兩部分組成,也即理想項(xiàng)和剩余項(xiàng):

理想項(xiàng)貢獻(xiàn)Cid p取自于實(shí)驗(yàn)關(guān)聯(lián)式,剩余項(xiàng)貢獻(xiàn)則是通過系綜平均獲得,其表達(dá)式如式(7)所示:

式中:KB為玻爾茲曼常數(shù),T和P分別為體系所設(shè)定的溫度和壓力,Hconf、Uconf分別為體系的構(gòu)型焓和構(gòu)型內(nèi)能(也即排除分子動(dòng)能和分子內(nèi)原子間的相互作用)。

3 結(jié)果與討論

3.1 密度(ρ)預(yù)測(cè)

采用4 種典型的分子力場(chǎng)模型(TraPPE、EPM2、Zhang、SAFT-γ Mie),采用分子動(dòng)力學(xué)模擬方法,在NPT 系綜下,對(duì)CO2流體在超臨界態(tài)工況條件下(T=600—900 K,p=30—100 MPa)的密度特性進(jìn)行了模擬預(yù)測(cè),且將模擬預(yù)測(cè)結(jié)果與NIST 數(shù)據(jù)庫(kù)中的參考數(shù)據(jù)進(jìn)行了對(duì)比,如圖1 所示。同時(shí),不同溫度下各分子力場(chǎng)密度預(yù)測(cè)的平均絕對(duì)相對(duì)誤差如表2所示。由圖1 和表2 可知,4 種力場(chǎng)模型對(duì)高溫高壓條件下CO2流體的密度預(yù)測(cè)值與NIST 參考值的平均絕對(duì)相對(duì)偏差均在3%以內(nèi),因此均具有較好預(yù)測(cè)精度。值得注意的是,3 種全原子力場(chǎng)(TraPPE、EPM2、Zhang)對(duì)密度的預(yù)測(cè)精度均隨著溫度的升高呈現(xiàn)逐漸減小的趨勢(shì),也即在高溫條件下具有更好的預(yù)測(cè)精度,而單點(diǎn)粗?；疭AFT-γ Mie 力場(chǎng)密度預(yù)測(cè)精度隨溫度的變化趨勢(shì)卻相反。盡管如此,簡(jiǎn)化的SAFT-γ Mie 力場(chǎng)總體上在超臨界態(tài)工況條件下(T=600—900 K,p=30—100 MPa)對(duì)密度具有相對(duì)最好的預(yù)測(cè)精度,這主要由于SAFT-γ Mie 力場(chǎng)模型參數(shù)的確定采用了一種全局化的參數(shù)優(yōu)化策略。

圖1 不同力場(chǎng)所預(yù)測(cè)的CO2 流體密度值與NIST 參數(shù)數(shù)據(jù)的對(duì)比Fig.1 Comparison of predicted CO2 fluid density values using different force fields and NIST parameter data

表2 不同力場(chǎng)在不同溫度下密度預(yù)測(cè)的平均絕對(duì)相對(duì)偏差Table 2 Average absolute relative deviation of density prediction using different force fields at different temperatures

3.2 比定壓熱容(Cp)預(yù)測(cè)

同樣在NPT 系綜下,對(duì)超臨界態(tài)工況下(T=600—900 K,p=30—100 MPa)CO2流體的比定壓熱容特性進(jìn)行了模擬預(yù)測(cè),模擬值與NIST 參考數(shù)據(jù)對(duì)比結(jié)果如圖2 所示。其中,表3 為不同溫度下各分子力場(chǎng)所預(yù)測(cè)的平均絕對(duì)相對(duì)誤差。由圖2 和表3 可知,針對(duì)超臨界態(tài)比定壓熱容,4 種力場(chǎng)模型所預(yù)測(cè)的平均絕對(duì)相對(duì)偏差均在1%以內(nèi),因此比密度的預(yù)測(cè)精度更高。由表3 可知,EPM2 和Zhang 力場(chǎng)模型對(duì)比定壓熱容的預(yù)測(cè)偏差均隨溫度的升高而逐漸減小,而TraPPE 和SAFT-γ Mie 力場(chǎng)的預(yù)測(cè)精度隨溫度的變化不顯著。總體上,在溫度T=600—900 K 和壓力p=30—100 MPa 的超臨界態(tài)范圍內(nèi),EPM2 力場(chǎng)對(duì)CO2流體的比定壓熱容特性具有最好的預(yù)測(cè)精度。

表3 超臨界工況下,不同力場(chǎng)在不同溫度下的比定壓熱容預(yù)測(cè)偏差Table 3 Average absolute relative deviation of specific heat capacity at constant pressure using different force fields at different temperatures

圖2 不同力場(chǎng)所預(yù)測(cè)的CO2 流體比定壓熱容與NIST 參數(shù)數(shù)據(jù)的對(duì)比Fig.2 Comparison of predicted CO2 specific heat capacity at constant pressure using different force fields and NIST parameter data

同時(shí),上述模擬結(jié)果表明,基于氣液相平衡特性擬合關(guān)聯(lián)的CO2力場(chǎng)模型對(duì)超臨界態(tài)熱力學(xué)性質(zhì)(如密度和比定壓熱容)同樣具有較好的預(yù)測(cè)精度。此外,盡管全原子力場(chǎng)模型具有更多的參數(shù)自由度,但在超臨界態(tài)熱力學(xué)性質(zhì)預(yù)測(cè)中,其預(yù)測(cè)性能相比粗?；?chǎng)模型并非具有顯著優(yōu)勢(shì)。由于在溫度T=600—900 K 和壓力p=30—100 MPa 的超臨界態(tài)區(qū)間內(nèi),CO2流體不具有顯著的變物性行為,因此4 種力場(chǎng)模型均體現(xiàn)出較好的預(yù)測(cè)精度。同時(shí),隨著溫度的升高,EPM2 和Zhang 兩種力場(chǎng)對(duì)密度和定壓比熱容的預(yù)測(cè)偏差均呈現(xiàn)出逐漸減小的趨勢(shì),這表明該力場(chǎng)的預(yù)測(cè)能力可推廣應(yīng)用于更高溫度的熱力工況。

4 結(jié)論

采用分子動(dòng)力學(xué)模擬方法結(jié)合4 種典型的半經(jīng)驗(yàn)型分子力場(chǎng),對(duì)超臨界態(tài)二氧化碳流體熱力學(xué)性質(zhì)進(jìn)行了模擬預(yù)測(cè)研究,并得出以下結(jié)論:

(1)盡管4 種力場(chǎng)的范德華作用參數(shù)通過擬合氣液相平衡數(shù)據(jù)而獲取,但在溫度T=600—900 K 和壓力p=30—100 MPa 的超臨界態(tài)區(qū)間內(nèi),4 種力場(chǎng)對(duì)流體密度和比定壓熱容特性均具有較好的預(yù)測(cè)精度,這表明4 種力場(chǎng)模型均具有較好的外延預(yù)測(cè)性能。

(2)EPM2 和Zhang 力場(chǎng)在對(duì)密度和比定壓熱容的模擬預(yù)測(cè)中,其預(yù)測(cè)偏差隨溫度升高均呈現(xiàn)出逐漸減下的趨勢(shì),這表明該力場(chǎng)可用于更高溫度工況條件下CO2流體熱力學(xué)性質(zhì)的預(yù)測(cè)。

(3)相比于簡(jiǎn)單且計(jì)算高效的單粒子SAFT-γ Mie 力場(chǎng),全原子力場(chǎng)模型盡管具有更多的參數(shù)自由度,但其在超臨界態(tài)區(qū)間內(nèi)對(duì)CO2流體熱力學(xué)特性的預(yù)測(cè)并不具有顯著優(yōu)勢(shì)。這表明全原子力場(chǎng)模型的力場(chǎng)參數(shù)對(duì)于超臨界態(tài)熱物性預(yù)測(cè)而言并非最優(yōu)解。