姜思雨, 郭 濤, 胡家文*

超臨界氮的四參數(shù)立方型狀態(tài)方程與分子動(dòng)力學(xué)模擬

姜思雨1, 2, 郭 濤1, 3, 胡家文1, 2*

(1. 河北地質(zhì)大學(xué) 河北省戰(zhàn)略性關(guān)鍵礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室, 河北 石家莊 050031; 2. 河北地質(zhì)大學(xué) 地球科學(xué)學(xué)院, 河北 石家莊 050031; 3. 河北地質(zhì)大學(xué)數(shù)理教學(xué)部, 河北 石家莊 050031)

根據(jù)精確的壓力–摩爾體積–溫度()數(shù)據(jù)及其外延結(jié)果, 本研究建立了超臨界氮(N2)的四參數(shù)立方型狀態(tài)方程。新方程最適宜的條件為273～4273 K、0～0.85 g·cm–3。在此范圍內(nèi), 新方程的最大體積偏差在0.41%以?xún)?nèi), 平均體積偏差在0.1%以?xún)?nèi); 等溫線的壓力范圍隨著溫度的上升而擴(kuò)大,最大壓力達(dá)2.9 GPa。當(dāng)密度增至1.03 g·cm–3、溫度增至5000 K時(shí), 體積偏差在0.6%～2.50%以?xún)?nèi), 最高壓力達(dá)3.3～5.0 GPa。新方程也可近似地應(yīng)用到密度為1.2 g·cm–3。此時(shí)相應(yīng)的溫壓可達(dá)5573 K、7.5 GPa。此外, 本研究還用分子動(dòng)力學(xué)方法模擬了超臨界氮在1473～5073 K、0.3～6.0 GPa范圍內(nèi)的性質(zhì)。所得結(jié)果以及現(xiàn)有的高溫或高密度實(shí)驗(yàn)數(shù)據(jù)、分子動(dòng)力學(xué)和蒙特卡洛(Monte Carlo)模擬數(shù)據(jù)(0.1～1.2 g·cm–3)均與新方程吻合得很好或較好。由新方程預(yù)測(cè)的逸度系數(shù)、剩余焓、剩余熵及其他熱力學(xué)性質(zhì)與參考模型的計(jì)算值均很吻合。這些結(jié)果表明新方程顯著地優(yōu)于現(xiàn)有的立方型、維里型和更復(fù)雜的地質(zhì)流體狀態(tài)方程。

超臨界氮; 狀態(tài)方程;性質(zhì); 逸度系數(shù); 剩余焓; 分子動(dòng)力學(xué)

0 引 言

N2不僅是太陽(yáng)系內(nèi)巨行星及太陽(yáng)系外行星的主要成分之一(Mochalov et al., 2010; Driver and Militzer, 2016; Dubois et al., 2017), 也是地球表面及內(nèi)部常見(jiàn)的流體組分之一。由于N2的臨界溫度(C)很低(126.192 K)、臨界壓力(C)也不太高(3.3958 MPa), 因此N2在自然界經(jīng)常以超臨界流體形式存在。大量研究結(jié)果表明, 含N2流體可以出現(xiàn)在不同的地質(zhì)環(huán)境中。例如, 含N2的流體包裹體常見(jiàn)于各種變質(zhì)巖(Fu et al., 2001; 張澤明等, 2006; Wang et al., 2018; 肖益林等, 2018)、幔源金剛石(Smith et al., 2014; Sobolev et al., 2019)、海相蒸發(fā)巖(Siemann and Ellendorff, 2001)、地?zé)崃黧w(Dubessy et al., 1989)、煤層(So?nicka and Lüders, 2019)和多種礦床(倪培等, 2018), 其中有的高密度包裹體在室溫下甚至呈現(xiàn)固態(tài)(Sobolev et al., 2019)。這些包裹體保存了地質(zhì)環(huán)境的許多信息, 通過(guò)對(duì)流體包裹體的組成、相變溫度(冰點(diǎn)、水合物的形成或熔化溫度、包裹體的部分或整體均一化溫度)等性質(zhì)的分析和相關(guān)的熱力學(xué)模擬可以獲取地質(zhì)過(guò)程的溫度、壓力和流體密度等重要物理化學(xué)條件及其變化值(盧煥章等, 2004), 進(jìn)而計(jì)算流體系統(tǒng)的氧逸度(O2)、pH值等(劉斌, 2011; Artimenko, 2017)。這些結(jié)果對(duì)于理解流體–礦物的相互作用、礦物、巖石和礦床(包括油氣藏)的形成與演化、變質(zhì)作用、地?zé)崤c巖漿活動(dòng)等許多過(guò)程十分重要(盧煥章等, 2004; 孫賀和肖益林, 2009; 倪培等, 2018; 肖益林等, 2018)。此外, 含N2的流體也會(huì)出現(xiàn)在煤層氣(Krooss et al., 1995; So?nicka and Lüders, 2019)、地幔捕虜體(Andersen et al., 1995; Berkesi et al., 2017)、熔體或巖漿(Yoshioka et al., 2018)、火山氣體與地?zé)崃黧w(Elkins et al., 2006; Fischer, 2008)。地質(zhì)流體中的N2可以以微量、少量或主量組分甚至純N2的形式出現(xiàn)(Andersen et al., 1993; Elvevold and Andersen, 1993; Krooss et al., 1995; Fu et al., 2001; 翟偉等, 2005; Marko et al., 2006; Ye et al., 2012; Smith et al., 2014; Yang et al., 2016; Wang et al., 2018)。這些流體及其地質(zhì)環(huán)境的溫壓條件可以從常溫、常壓變化到1300 ℃以上、6～8 GPa (Van den Kerkhof et al., 1991; Fu et al., 2001; Xiao et al., 2002; 張澤明等, 2006; Smith et al., 2014; Sobolev et al., 2019)。

在地球或天體系統(tǒng)內(nèi), 分子態(tài)N2的主要存在形式是流體。此外, 氮還有許多其他存在形式, 包括無(wú)機(jī)、有機(jī)的離子或分子形式(固態(tài)或流體態(tài))。其中, 分子態(tài)的流體氮與氮的其他形式之間通過(guò)遷移和反應(yīng)相互轉(zhuǎn)化(Roskosz et al., 2013; Johnson and Goldblatt, 2015; Sokol et al., 2017; Mallik et al., 2018)。因此, 要正確地理解含氮物種之間的轉(zhuǎn)化(包括相平衡和化學(xué)平衡)及深部氮循環(huán), 必須掌握含氮物質(zhì)的熱力學(xué)性質(zhì)(Roskosz et al., 2013; Yang et al., 2016; Mikhail et al., 2017; Sokol et al., 2017; Mallik et al., 2018; Chen et al., 2019)。

N2也是天然氣的主要成分之一。天然氣工業(yè)的許多環(huán)節(jié)都需要了解N2及相關(guān)混合物的熱力學(xué)性質(zhì)。N2是許多化工生產(chǎn)、超臨界水氧化、廢水處理、生物質(zhì)加工、油田后期開(kāi)采過(guò)程的原料或產(chǎn)物(周健等, 1999; 楊光璐等, 2004; García-Jarana et al., 2013; Kraussler et al., 2016), 也是許多含氮物質(zhì)燃燒、爆炸或沖擊壓縮過(guò)程的主要產(chǎn)物(Driver and Militzer, 2016; Dubois et al., 2017; Fu et al., 2019)。

對(duì)上述各種過(guò)程的理解、模擬、設(shè)計(jì)和控制需要準(zhǔn)確地掌握N2及其混合物的壓力–摩爾體積–溫度–組成()、熱容、焓、熵及許多其他的熱力學(xué)性質(zhì), 其中和熱容是最基本的性質(zhì)。和熱容等基礎(chǔ)數(shù)據(jù)被總結(jié)成各種熱力學(xué)模型, 如狀態(tài)方程、Helmholtz和Gibbs自由能模型。這些模型與微積分運(yùn)算相結(jié)合可以推算出其他熱力學(xué)性質(zhì)。

1 N2的PVT實(shí)驗(yàn)或模擬數(shù)據(jù)

迄今為止, 人們已通過(guò)實(shí)驗(yàn)手段獲得流體N2的大量熱力學(xué)數(shù)據(jù), 包括幾組極高溫壓下沖擊壓縮實(shí)驗(yàn)的數(shù)據(jù)(Radousky and Ross, 1988; Nellis et al., 1991; Mochalov et al., 2010), 而其余實(shí)驗(yàn)數(shù)據(jù)的最高溫壓只有1800 K和2.2 GPa。Span et al. (2000)對(duì)截止到1997年的N2的各種熱力學(xué)實(shí)驗(yàn)數(shù)據(jù)做了比較全面的搜集、統(tǒng)計(jì)和不確定度評(píng)估, 盡管也有一些實(shí)驗(yàn)數(shù)據(jù)被他們遺漏(Malbrunot and Vodar, 1973; Antanovich and Plotnikov, 1977; Roder, 1981; Belak et al., 1988; Perkins et al., 1991; 黃榮榮, 1995)。在近20年中, 陸續(xù)有一些N2的實(shí)驗(yàn)數(shù)據(jù)發(fā)表, 但均局限在中低溫壓范圍(Mourato et al., 1999; Evers et al., 2002; Patil, 2005; Saleh and Wendland, 2005; Chamorro et al., 2006; Atilhan, 2007; Mantilla et al., 2010; Sakoda et al., 2012; Yang et al., 2015; 徐航等, 2015; Cheng et al., 2019)。總體而言, 絕大多數(shù)流體N2的熱力學(xué)實(shí)驗(yàn)數(shù)據(jù)都在1000 K、1 GPa以下。其中相當(dāng)一部分?jǐn)?shù)據(jù)存在較大的不確定度或系統(tǒng)偏差; 而且,隨著溫度或壓力的升高, 不確定度或系統(tǒng)偏差增加較快。對(duì)于許多高溫高壓生產(chǎn)或相關(guān)研究(如燃燒、含能材料、爆轟與沖擊波物理、地球與行星深部的物理化學(xué)研究)而言, 這種情況遠(yuǎn)不能滿(mǎn)足實(shí)際需要。

為了獲取更多的高溫高壓熱力學(xué)數(shù)據(jù), 人們采用原子–分子尺度上的計(jì)算機(jī)模擬, 如分子動(dòng)力學(xué)(molecular dynamics, MD)、蒙特卡洛(monte carlo, MC)、密度泛函理論(density function theory, DFT)、第一性原理(first principles, FP)模擬、量子力學(xué)從頭計(jì)算法(), 或這些方法的某種聯(lián)合應(yīng)用。這些方法只需少量參數(shù)或者無(wú)需任何經(jīng)驗(yàn)參數(shù)就可以獲得很寬溫壓范圍內(nèi)(包括極端溫壓下)流體的單相、多相系統(tǒng)和化學(xué)反應(yīng)的各種熱力學(xué)數(shù)據(jù), 從而極大地拓寬了流體熱力學(xué)數(shù)據(jù)的溫壓范圍, 同時(shí)也顯示了很好的預(yù)測(cè)能力(Johnson et al., 1984; Etters et al., 1986; Kress et al., 2000; Mazevet et al., 2001; Brennan and Rice, 2002; Krukowski and Str?k, 2006; Maiti et al., 2007; Str?k and Krukowski, 2007; Boates and Bonev, 2009; Driver and Militzer, 2016; Lv et al., 2018; Fu et al., 2019)。與傳統(tǒng)的高溫、高壓實(shí)驗(yàn)方法相比, 這類(lèi)方法在時(shí)間、成本、安全性和適用范圍等方面具有巨大的優(yōu)勢(shì)。

目前, 通過(guò)實(shí)驗(yàn)(Radousky and Ross, 1988; Nellis et al., 1991; Mochalov et al., 2010)或模擬(Johnson et al., 1984; Etters et al., 1986; Kress et al., 2000; Mazevet et al., 2001; Brennan and Rice, 2002; Maiti et al., 2007; Str?k and Krukowski, 2007; Boates and Bonev, 2009; Driver and Militzer, 2016; Fu et al., 2019)方法已獲得了不少高溫、高壓(包括極端溫壓)條件下N2的熱力學(xué)數(shù)據(jù)。對(duì)于地球與行星深部研究而言, 1000～5000 K、1～5 GPa范圍內(nèi)流體N2的熱力學(xué)數(shù)據(jù)特別重要, 但相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)卻很少。相關(guān)的模擬數(shù)據(jù)雖然較多一些, 但有不少數(shù)據(jù)具有顯著的系統(tǒng)偏差(Johnson et al., 1984; Driver and Militzer, 2016), 而且很多模擬數(shù)據(jù)只有圖形報(bào)道(Etters et al., 1986; Str?k and Krukowski, 2007), 難以精確地再現(xiàn)(Maiti et al., 2007)或報(bào)道有誤(Krukowski and Str?k, 2006)。因此, 目前仍然需要補(bǔ)充高溫、高壓下N2的熱力學(xué)數(shù)據(jù)。

2 N2的熱力學(xué)模型

自20世紀(jì)70年代以來(lái), 許多研究致力于發(fā)展地質(zhì)流體(包括含N2流體)的熱力學(xué)模型, 其中主要是狀態(tài)方程(Bakker, 2003; Gottschalk, 2007; 段振豪, 2010; Zhang et al., 2016)。目前, 含N2流體的熱力學(xué)模型大致可以分為4類(lèi):

①立方型狀態(tài)方程。自從van der Waals (vdW)提出第一個(gè)立方型狀態(tài)方程以來(lái), 不斷有新的立方型狀態(tài)方程發(fā)表。Redlich and Kwong (1949)(RK49)將vdW方程的引力參數(shù)改為/1/2, 所得方程如下:

式中:為協(xié)體積(co-volume);為通用氣體常數(shù);為摩爾體積。此方程在中低壓(幾十MPa)下的精度得到了顯著的提高, 但在中高壓下仍然欠佳。盡管如此, Holloway (1977)直接將該方程用于除H2O和CO2以外的非極性或弱極性地質(zhì)流體組分(包括N2)。為了計(jì)算含CH4和N2的地質(zhì)流體的熱力學(xué)性質(zhì), Bakker (1999)(B99)將RK49方程的參數(shù)改為溫度的線性函數(shù):

式中:0和1為常量;0= 273.15 K。此方程中和單位分別為K、MPa和cm3·mol–1。對(duì)于N2,= 27 cm3·mol?1。Bakker將此方程應(yīng)用到1300 K、1 GPa, 但偏差較大。

總體而言, 這類(lèi)方程形式簡(jiǎn)單, 參數(shù)較少, 在給定溫壓之下有解析形式的體積根, 一般都具有較好的預(yù)測(cè)性, 便于以同樣的形式推廣到多組分混合物, 而且所需要的混合參數(shù)較少, 應(yīng)用方便。目前, 這類(lèi)方程數(shù)以百計(jì), 在各種工業(yè)流體的熱力學(xué)計(jì)算或模擬中得到了廣泛的應(yīng)用。不過(guò), 立方型方程所適用的溫壓范圍一般都在幾百K和100 MPa以?xún)?nèi), 在中高密度條件下對(duì)摩爾體積及相關(guān)廣延性質(zhì)的計(jì)算偏差通常較大或很大。

②維里型狀態(tài)方程及其衍生形式。此類(lèi)方程也較多, 但其中很少有方程能在很寬的–范圍內(nèi)保持較高的精度。例如, Lee and Kesler (1975)和Soave (1995, 1999)方程主要面向工業(yè)流體, 不適合高壓地質(zhì)流體。Saxena and Fei (1987a, 1987b)提出兩種對(duì)應(yīng)態(tài)形式的維里型方程, 其中前者Saxena and Fei (1987a)采用了分段函數(shù), 后者Saxena and Fei (1987b)在中低壓下不正確。但它們的偏差都較大或很大。Holland and Powell (1991)(HP91)將RK49方程近似地轉(zhuǎn)化為以和為自變量的分式方程, 并加上2個(gè)補(bǔ)償項(xiàng):

式中:、和的單位分別為K、kbar、kJ·kbar?1·mol?1(=10 cm3·mol?1)和kJ2·kbar?1·K1/2·mol?2,= 0.0083143 kJ·K?1·mol?1。不過(guò), 他們并未報(bào)道其方程對(duì)N2的具體效果。

Duan et al. (1992)(DMW92)利用CH4的實(shí)驗(yàn)數(shù)據(jù)和MD模擬結(jié)果, 針對(duì)非極性的地質(zhì)流體組分(包括N2)建立了一種對(duì)應(yīng)態(tài)形式的通用維里型狀態(tài)方程:

Belonoshko and Saxena (1992)(BS92)利用MD模擬結(jié)果和exp-6勢(shì)能參數(shù)建立了一種對(duì)應(yīng)態(tài)形式的維里型狀態(tài)方程:

③半經(jīng)驗(yàn)的統(tǒng)計(jì)力學(xué)狀態(tài)方程。此類(lèi)方程中有一些不太復(fù)雜的方程(Christoforakos and Franck, 1986; Heilig and Franck, 1989), 但多數(shù)都比較復(fù)雜或很復(fù)雜, 如目前被廣泛應(yīng)用的有統(tǒng)計(jì)締合流體理論(statistical associating fluid theory, SAFT)、微擾理論(perturbation theory, PT)狀態(tài)方程(Churakov and Gottschalk, 2003; Sun et al., 2015), 其中SAFT方程一般所需參數(shù)較少, 預(yù)測(cè)性較好, 但對(duì)密度或體積的計(jì)算偏差較大。這類(lèi)方程主要用于工業(yè)流體, 相應(yīng)的溫壓范圍比較有限, 一般不適于高溫、高壓下的地質(zhì)流體。

④Helmholtz自由能模型。其中有一類(lèi)是比較復(fù)雜的專(zhuān)用模型, 如Span et al. (2000)和Jacobsen et al. (1986)的Helmholtz自由能模型。這類(lèi)模型適用溫壓范圍較寬, 也很精確, 但需要擬合的參數(shù)很多, 很難以同樣的形式推廣到混合物。另一類(lèi)是經(jīng)過(guò)簡(jiǎn)化、具有一定通用性的Helmholtz自由能模型, 例如Sun and Ely (2004)與Mao et al. (2017)的模型。這類(lèi)模型比較精確, 但比維里型方程要稍復(fù)雜一些, 能夠應(yīng)用的溫壓范圍比較有限, 不適合于高溫、高壓地質(zhì)流體。

關(guān)于含N2地質(zhì)流體的上述熱力學(xué)模型一般均有較大的偏差, 而且實(shí)際可用的溫壓范圍與原作者所聲明的結(jié)果經(jīng)常有很大(或較大)的差別。這也印證了本研究組近期關(guān)于地質(zhì)流體狀態(tài)方程的研究結(jié)果(王艷哲, 2015; 郭濤等, 2016; 郭磊, 2017)。

目前, 關(guān)于含N2地質(zhì)流體的熱力學(xué)模型雖然很多, 但簡(jiǎn)單、精確、適用范圍寬廣的熱力學(xué)模型仍然比較缺乏。鑒于目前地質(zhì)流體組分及混合物的高溫、高壓熱力學(xué)數(shù)據(jù)仍然短缺, 采用比較簡(jiǎn)單的模型來(lái)模擬地質(zhì)流體的熱力學(xué)性質(zhì)是一種明智的選擇。

3 四參數(shù)立方型狀態(tài)方程的建立

考慮到立方型狀態(tài)方程具有多方面的優(yōu)點(diǎn), 在適用范圍的擴(kuò)展方面仍有很大的潛力, 本研究擬建立一種新的立方型狀態(tài)方程, 力圖在較寬的溫壓范圍內(nèi)精確地再現(xiàn)或預(yù)測(cè)含N2流體的多種熱力學(xué)性質(zhì)。為此, 首先需要建立純N2流體的狀態(tài)方程。

3.1 方程的表達(dá)式

為了兼顧簡(jiǎn)潔性、精確度和適用范圍, 本文采用以下的四參數(shù)立方型狀態(tài)方程:

式中:2、3和均為溫度的函數(shù);為常數(shù);= 8.31451 J·K–1·mol–1;、與的單位分別是MPa、cm3·mol?1和K。在Lennard-Jones流體第二維里系數(shù)的統(tǒng)計(jì)力學(xué)解析表達(dá)式中, 溫度的指數(shù)均為1/4的整數(shù)倍。據(jù)此, 胡家文(2002)采用溫度指數(shù)為1/4的整數(shù)倍(或相近值)的簡(jiǎn)單多項(xiàng)式來(lái)擬合第二維里系數(shù), 效果很好。郭濤等(2016)關(guān)于水的立方型狀態(tài)方程中引力參數(shù)采用了類(lèi)似的形式。結(jié)果表明, 這種表達(dá)式可以有效地外延到極高的溫度。

對(duì)方程(17)中2、3與的關(guān)系, 本研究嘗試了不同的表達(dá)式(包括改變表達(dá)式的形式或增減參數(shù)的數(shù)量), 結(jié)果發(fā)現(xiàn)以下的形式效果最好:

3.2 方程的參數(shù)擬合

Span et al. (2000)從截止到1997年的N2的各種熱力學(xué)數(shù)據(jù)中選用最好的數(shù)據(jù)建立了一種高精度的Helmholtz自由能模型。該模型最適宜的溫壓應(yīng)用范圍是63.151～1000 K、0～2.2 GPa。從三相點(diǎn)到523 K、0～12 MPa或240～523 K、0～30 MPa, 密度的不確定度為0.02%; 當(dāng)壓力很高(>1 GPa)時(shí), 密度的不確定度可能會(huì)達(dá)到0.6%; 在270～350 K、12 MPa以下, 密度的不確定度僅為0.01%。鑒于該模型的精度和可靠性, 美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究所的官方網(wǎng)站將它推薦為N2的參考模型。不過(guò), 該模型的形式很復(fù)雜, 線性和非線性參數(shù)都很多, 不便于擬合或應(yīng)用, 也難以同樣的形式應(yīng)用于混合物, 因而其主要作用在于為科學(xué)研究提供可靠的數(shù)據(jù), 有時(shí)也被用作混合物熱力學(xué)模型中的純組分模型。

為了保證新方程的準(zhǔn)確性及外延結(jié)果的可靠性, 方程(17)參數(shù)擬合時(shí)所用的數(shù)據(jù)全部來(lái)自上述參考模型。為了使方程(17)具有盡可能寬的適用范圍, 本研究對(duì)參考模型在上述范圍內(nèi)的計(jì)算結(jié)果向高溫、高壓方向進(jìn)行了外延。

對(duì)于大多數(shù)流體而言, 中低密度時(shí)的等容線非常接近于直線; 在中高密度條件下, 等容線在中低溫段會(huì)有明顯的彎曲, 但在高溫段彎曲程度會(huì)明顯減弱。據(jù)此很容易對(duì)等容線進(jìn)行精確的擬合, 同時(shí)也適于做大幅度的外延。本研究對(duì)N2等容線的擬合與外延均采用如下的多項(xiàng)式:

外延結(jié)果與參考模型在2000 K、2.2 GPa以下的計(jì)算結(jié)果一起用于方程(17)的參數(shù)擬合。

由于流體在臨界點(diǎn)附近具有奇異性, 臨界點(diǎn)附近性質(zhì)的擬合難度遠(yuǎn)高于遠(yuǎn)離臨界點(diǎn)的超臨界區(qū)。這是因?yàn)閷?duì)臨界奇異性的嚴(yán)格處理需要采用非解析的(即非經(jīng)典的)理論模型或方法, 如重整化群和跨接函數(shù)法, 但這些模型或方法極為復(fù)雜, 計(jì)算量也很大(郭濤等, 2016)。雖然一些半經(jīng)驗(yàn)的方法可以近似地表達(dá)流體的臨界奇異性, 但它們主要針對(duì)中低溫壓下的工業(yè)流體。由于上述原因, 加之N2的臨界溫度很低, 新方程的參數(shù)擬合與優(yōu)化未考慮臨界點(diǎn)附近的區(qū)域。

經(jīng)過(guò)反復(fù)嘗試, 新方程的參數(shù)擬合與優(yōu)化能夠確保方程精確地適用于273.15～4273.15 K?？紤]到一些高溫、高壓研究的需要, 參數(shù)的優(yōu)化力求讓新方程能在高達(dá)5000 K左右仍然具有較高的精度。為了兼顧流體包裹體研究的需要, 參數(shù)的優(yōu)化力求使所得方程能較好地外延到200 K附近。這一溫度比CO2的三相點(diǎn)(216.592 K)還要低十幾度。

新方程的參數(shù)擬合采用非線性最小二乘法。經(jīng)過(guò)反復(fù)優(yōu)化, 所得結(jié)果見(jiàn)表1。

注: E+表示×10+n, E?表示×10?n。

3.3 對(duì)參考模型及其高溫外延結(jié)果的再現(xiàn)

表1中的參數(shù)可以精確地再現(xiàn)參考模型的計(jì)算結(jié)果及其等容線的高溫外延結(jié)果。根據(jù)再現(xiàn)結(jié)果, 方程(17)最適宜的條件是273.15～4273.15 K、0～0.85 g·cm?3。方程(17)在273.15～2000 K區(qū)間內(nèi)的再現(xiàn)結(jié)果見(jiàn)圖1, 在2000～4273.15 K區(qū)間內(nèi)的高溫外延的再現(xiàn)結(jié)果見(jiàn)圖2。

圖1中的最大體積偏差(the maximum deviation of molar volume, MDV)為0.306%, 平均體積偏差(average deviation of molar volume, ADV)為0.089%?？傮w上看, 在450 K以下, 隨著溫度的下降, 體積偏差增加較快; 在450 K以上, 隨著溫度的上升, 體積偏差緩慢下降。圖2中, 方程(17)的MDV=0.402%, ADV=0.091%。

cal代表方程(17); ref代表Span et al. (2000)的參考模型。

圖a: 1.=33 cm3·mol?1; 2.=34 cm3·mol?1; 3.=35 cm3·mol?1; 4.=38 cm3·mol?1; 5.=40 cm3·mol?1; 6.=45 cm3·mol?1; 7.=50 cm3·mol?1; 8.=55 cm3·mol?1; 9.=60 cm3·mol?1; 10.=80 cm3·mol?1; 11.=100 cm3·mol?1; 圖b: 1.=200 cm3·mol?1; 2.=500 cm3·mol?1; 3.= 1000 cm3·mol?1; 4.=2000 cm3·mol?1; 5.=4000 cm3·mol?1; 參考模型數(shù)值據(jù)Span et al., 2000。

圖2 方程(17)與參考模型等容外延結(jié)果的比較(=2000～4273.15 K,0～0.85 g·cm?3)

Fig.2 Comparison of Equation (17) and the results of isochoric extrapolation of the reference model

3.4 方程的適用范圍

參考模型在4273.15 K以上的外延結(jié)果也與方程(17)吻合得很好(圖3a), 其中MDV=0.60%, ADV= 0.27%。方程(17)向低溫方向適當(dāng)?shù)耐庋咏Y(jié)果如圖3(b)所示, 其中MDV<2%。

由于立方型方程能夠應(yīng)用的最高密度max比較有限, 當(dāng)密度接近或超過(guò)max時(shí), 方程的偏差會(huì)隨著密度的增大而迅速增大。當(dāng)密度范圍一定時(shí), 方程的最大偏差(或偏差的變化范圍)也是一定的。根據(jù)方程(17)的最大偏差, 本研究劃分出幾個(gè)密度區(qū)間(表2), 并按區(qū)間統(tǒng)計(jì)ADV與MDV, 其中也包括273.15 K以下和4273.15 K以上的偏差統(tǒng)計(jì)結(jié)果。

在273.15～4273.15 K區(qū)間內(nèi), 當(dāng)max為0.850、0.926、0.964、1.000和1.030 g·cm?3時(shí), 方程(17)相對(duì)于參考模型的MDV分別為0.41%、0.95%、1.39%、2.00%和2.50%。為了便于應(yīng)用, 本文用方程(17)計(jì)算了這5個(gè)密度對(duì)應(yīng)的等容線(圖4), 并用經(jīng)驗(yàn)式(21)對(duì)它們進(jìn)行了擬合(表3)。

對(duì)于每條等容線, 壓力隨著溫度的增大而增大。因此, 在方程能夠應(yīng)用的max以?xún)?nèi), 溫度越低, 方程能夠應(yīng)用的最大壓力(max)也就越低。盡管如此, 方程(17)的適用范圍基本上涵蓋了一些典型的地溫、地壓梯度線(Tarney and Windley, 1977; Lambert, 1983; Chapman, 1986)及典型的變質(zhì)巖相的穩(wěn)定區(qū)(Dipietro, 2013), 包括高溫、超高溫變質(zhì)作用的溫壓條件(劉守偈等, 2008; 魏春景, 2016; 李旭平等, 2019)。因此, 對(duì)于不同的地殼溫壓梯度, 方程(17)可適用于下地殼乃至上地幔的溫壓條件。

4 新方程的檢驗(yàn)

4.1 與高溫、高壓實(shí)驗(yàn)數(shù)據(jù)的比較

現(xiàn)有的高溫、高壓實(shí)驗(yàn)數(shù)據(jù)并不多。本文選取了5組具有代表性的高溫、高壓實(shí)驗(yàn)數(shù)據(jù)(Tsiklis and Polyakov, 1968; Robertson and Babb, 1969; Malbrunot and Vodar, 1973; Mills et al., 1975; Antanovich and Plotnikov, 1976)與新方程進(jìn)行了比較, 結(jié)果如圖5和表4所示。其中Robertson and Babb (1969)的數(shù)據(jù)中約有1/3超過(guò)了方程(17)最適宜的密度范圍(0～0.85 g·cm?3); Mills et al. (1975)的數(shù)據(jù)大部分都超過(guò)了0.85 g·cm?3, 所以在此只能選取密度在0.85 g·cm?3以下及附近(0.80～0.97 g·cm?3)的數(shù)據(jù)。這些數(shù)據(jù)相對(duì)于參考模型計(jì)算值的偏差見(jiàn)表4。其中Malbrunot and Vodar (1973)的體積數(shù)據(jù)在700 K以上有明顯的負(fù)偏差, 并且隨著溫度的升高而增大, 最大值達(dá)?5%, 平均約為?2%。

從總體上看, 新方程在其適用范圍內(nèi)與實(shí)驗(yàn)數(shù)據(jù)均吻合得很好。其中Malbrunot and Vodar (1973)在800 K以上的體積數(shù)據(jù)由于有負(fù)偏差而均明顯地低于方程(17)的計(jì)算值。與Mills et al. (1975)的數(shù)據(jù)相比, 方程(17)的偏差總體上較大。這是因?yàn)槠涠鄶?shù)數(shù)據(jù)的壓力(或密度)均顯著地超過(guò)了方程(17)最適宜的范圍。不過(guò), 隨著溫度的升高或壓力的降低, 方程(17)的偏差下降得很明顯。這與方程(17)對(duì)擬合數(shù)據(jù)的再現(xiàn)偏差的變化趨勢(shì)是一致的。

cal代表方程(17); ref代表參考模型; 參考模型數(shù)據(jù)據(jù)Span et al., 2000。

表2 方程(17)的體積偏差與溫度、密度范圍的關(guān)系

注: ADV為平均體積偏差, MDV為最大體積偏差。

4.2 與高溫、高壓MD、MC模擬數(shù)據(jù)的比較

目前, 雖然原子–分子尺度上的計(jì)算機(jī)模擬提供了很多高溫、高壓下N2的數(shù)據(jù), 但其中只有少數(shù)幾組數(shù)據(jù)可與新方程做比較。例如, Etters et al. (1986)將的數(shù)值計(jì)算結(jié)果擬合為函數(shù), 同時(shí)加上多極相互作用和分子內(nèi)相互作用的貢獻(xiàn), 得到了N2的分子間勢(shì)能函數(shù), 并根據(jù)這一函數(shù)用恒壓MC方法模擬了N2的性質(zhì)(300～1500 K, 0.2～15 GPa)。Str?k and Krukowski (2007)采用高精度的量子力學(xué)從頭算方法得到了N2的分子間勢(shì)能, 并將這些勢(shì)能用于N2的性質(zhì)的MD模擬(300～2000 K, 0～40 GPa)。這兩組數(shù)據(jù)都有很多高密度數(shù)據(jù), 但本文只需要其中≤1.2 g·cm?3的數(shù)據(jù), 其范圍見(jiàn)表5。

Etters et al. (1986)數(shù)據(jù)中0.85 g·cm?3的數(shù)據(jù)點(diǎn)約占2/3, 其中最大密度達(dá)1.17 g·cm?3; Str?k and Krukowski (2007)的數(shù)據(jù)中0.85 g·cm?3的數(shù)據(jù)點(diǎn)也不少, 其中最大密度達(dá)到了1.20 g·cm?3。

上述兩組模擬數(shù)據(jù)的最高壓力為4.43 GPa, 而最高溫度只有2000 K。為了檢驗(yàn)方程(17)在高溫、高壓下的精度, 本研究應(yīng)用Materials Studio軟件的Forcite模塊在更寬的-范圍內(nèi)(1473.15～5073.15 K, 0.3～6.0 GPa)對(duì)中高密度(0.4～1.196 g·cm?3)的N2進(jìn)行了MD模擬。

表3 等容線方程(21)的參數(shù)

注: E+表示×10+n。

(a) 數(shù)據(jù)據(jù)Tsiklis and Polyakov, 1968; (b) 數(shù)據(jù)據(jù)Robertson and Babb, 1969; (c) 數(shù)據(jù)據(jù)Malbrunot and Vodar, 1973; (d) 數(shù)據(jù)據(jù)Antanovich and Plotnikov, 1976; (e) 數(shù)據(jù)據(jù)Mills et al., 1975。

模擬采用系綜(固定分子數(shù)的等溫等壓系綜), 其對(duì)象為3000個(gè)N2分子, 其中N2的分子間勢(shì)能取自COMPASS 庫(kù)(Yang et al., 2000)。每個(gè)宏觀狀態(tài)點(diǎn)的模擬分成幾個(gè)階段, 微觀狀態(tài)之間的時(shí)間步長(zhǎng)均為1 fs。第一階段是初步模擬, 時(shí)長(zhǎng)為50 ps。實(shí)際上系統(tǒng)趨于熱力學(xué)平衡只需要20～40 ps。在后面的每個(gè)階段, 時(shí)長(zhǎng)都是1000 ps, 共計(jì)106個(gè)微觀狀態(tài), 其中初態(tài)為前一階段模擬的終態(tài)。每個(gè)階段的模擬結(jié)束后,和的系綜平均值及不確定度都會(huì)被記錄。如果最近兩個(gè)階段的平均值之差在前一階段的不確定度范圍內(nèi), 則停止模擬。最后結(jié)果取最后一個(gè)階段的平均值。在中高密度(>0.3 g?cm?3)條件下, 前三個(gè)階段的計(jì)算即足以達(dá)到模擬要求。模擬的范圍見(jiàn)表5。比較結(jié)果表明, 這些數(shù)據(jù)與參考模型或其高溫外延結(jié)果均吻合得很好。

從圖6可以看出, 表5中的3組模擬結(jié)果與新方程均吻合得較好。其中較大的偏差都出現(xiàn)在高密度區(qū), 相應(yīng)的ADV和MDV也列于表5。

表4 一些高溫、高壓實(shí)驗(yàn)數(shù)據(jù)及其相對(duì)于參考模型的體積偏差(參考模型據(jù)Span et al., 2000)

注: ADV為平均體積偏差; MDV為最大體積偏差;為數(shù)據(jù)點(diǎn)的數(shù)量。

表5 方程(17)與一些高溫、高壓PVT或PρT模擬數(shù)據(jù)的比較

注: ADV為平均體積偏差; MDV為最大體積偏差;為數(shù)據(jù)點(diǎn)的數(shù)量。

圖6 方程(17)與分子動(dòng)力學(xué)、Monte Carlo模擬數(shù)據(jù)的比較(圖a數(shù)據(jù)據(jù)Etters et al., 1986; 圖b數(shù)據(jù)據(jù)Str?k and Krukowski, 2007)

5 與其他地質(zhì)流體狀態(tài)方程的比較

5.1 與參考模型及其高溫外延結(jié)果的比較

HP91方程并沒(méi)有說(shuō)明其適用范圍。此方程與參考模型及其高溫外延結(jié)果的比較見(jiàn)圖7。結(jié)果表明, 影響HP91方程精度的主要因素是溫度, 這是因?yàn)閰?shù)與的關(guān)系過(guò)于簡(jiǎn)單。在400 K以上, 該方程的體積偏差隨著溫度的上升以近似線性的方式快速增大, 到1300 K附近達(dá)4%左右(圖7a)。這些結(jié)果表明, HP91方程的應(yīng)用不宜超過(guò)1300 K。在1300 K以下, 該方程與參考模型的比較見(jiàn)圖7b。其中, 體積偏差隨著壓力的增加有明顯的下降趨勢(shì)。

據(jù)報(bào)道, DMW92方程可用到2000 K、2.0 GPa, 平均偏差約為1.5%。為了厘清DMW92方程實(shí)際的適用范圍, 本研究將該方程與參考模型及其高溫外延結(jié)果進(jìn)行了比較(圖8)。結(jié)果表明: ①壓力是影響該方程精確度的主要因素, 其中壓力與體積偏差之間具有顯著的線性相關(guān)性; ②當(dāng)壓力增至0.7 GPa時(shí), MDV約為3%左右, 而相應(yīng)的溫度可達(dá)4273.15 K左右。當(dāng)壓力增至1.0 GPa時(shí), MDV增大到5%左右。

BS92方程據(jù)稱(chēng)可用于400～4000 K、0.5～100 GPa。該方程與參考模型及其高溫外延結(jié)果的比較如圖9所示。很明顯, 在0.8 GPa以下, 體積偏差隨著壓力的下降迅速增大, 其中MDV達(dá)9.54%; 在1.2 GPa以上, 體積偏差隨著壓力的增大也上升較快, 在2.9 GPa附近達(dá)到6.40%。

據(jù)原文報(bào)道, B99方程可應(yīng)用到1300 K和1 GPa。本文用參考模型及其高溫外延結(jié)果對(duì)它進(jìn)行了驗(yàn)證。圖10可見(jiàn), 體積偏差隨著壓力的上升以近似于直線的方式快速增大, 在0.75 GPa和1 GPa時(shí)達(dá)到5.01%和7.34%。此方程可以應(yīng)用的溫度遠(yuǎn)高于1300 K,但必須限制壓力的范圍。

5.2 與高溫、高壓MD、MC模擬數(shù)據(jù)的比較

HP91、DMW92、BS92與B99方程與表5中的MD、MC模擬數(shù)據(jù)的比較結(jié)果見(jiàn)圖11。比較時(shí)對(duì)各方程及模擬數(shù)據(jù)的限制條件如下, HP91:≤1300 K,≤2.5 GPa; BS92:≥0.5 GPa; B99:≤1.0 GPa; DMW92:≤2.0 GPa; 方程(17)沒(méi)有限制。這些條件要么在各方程原文聲明的范圍內(nèi), 要么在本研究所確認(rèn)的適用范圍內(nèi)。

圖7 HP91方程與參考模型的對(duì)比(cal數(shù)據(jù)據(jù)Holland and Powell, 1991; ref數(shù)據(jù)據(jù)Span et al., 2000)

圖8 DMW92方程與參考模型及其高溫外延結(jié)果的比較(cal數(shù)據(jù)據(jù)Duan et al., 1992; ref數(shù)據(jù)據(jù)Span et al., 2000)

圖9 BS92方程與參考模型及其高溫外延結(jié)果的比較(cal數(shù)據(jù)據(jù)Belonoshko and Saxena, 1992; ref數(shù)據(jù)據(jù)Span et al., 2000)

圖10 B99方程與參考模型及其高溫外延結(jié)果的比較(cal數(shù)據(jù)據(jù)Bakker, 1999; ref數(shù)據(jù)據(jù)Span et al., 2000)

從總體上看, HP91、DMW92與B99方程的偏差均較大或很大, 方程(17)的偏差最小, 其次是BS92方程。與Etters et al. (1986)、Str?k and Krukowski (2007)和本研究的模擬結(jié)果相比, 方程(17)的(ADV, MDV)分別為(1.40%, 3.89%)、(0.81%, 2.99%)和(1.18%, 3.52%); BS92方程的(ADV, MDV)分別為(2.17%, 4.70%)、(1.08%, 3.82%)和(1.99%, 5.02%), 都明顯大于方程(17)的相應(yīng)偏差, 其主要原因是BS92方程在1 GPa以下的偏差較大。

圖11進(jìn)一步印證了方程(17)及上述4個(gè)方程與參考模型及其高溫外延結(jié)果的比較結(jié)果。這也表明, 方程(17)可以擴(kuò)展到比圖4中的等容線更高的條件, 但必然伴隨著偏差的上升。

根據(jù)方程(17)計(jì)算的等容線, 當(dāng)≤1.2 g·cm–3時(shí), 方程(17)在200～5573.15 K范圍內(nèi)的結(jié)果都是合理的。結(jié)合本研究的結(jié)果可以推測(cè), 當(dāng)≤1.2 g·cm–3時(shí), 如果≤5073.15 K, 方程(17)的MDV約為4%, 相應(yīng)的壓力不超過(guò)6.8 GPa; 如果≤5573.15 K, 方程(17)的MDV約為4.5%, 相應(yīng)的壓力不超過(guò)7.5 GPa。

以上結(jié)果表明, 方程(17)在其適用范圍內(nèi)均明顯優(yōu)于HP91、DMW92、BS92與B99方程。其中DMW92與B99方程的適用范圍與原文報(bào)道的結(jié)果相差很大, 而B(niǎo)S92方程在0.8 GPa以下偏差較大。

6 其他熱力學(xué)性質(zhì)的預(yù)測(cè)

含氮物質(zhì)的遷移、交換、轉(zhuǎn)化過(guò)程(如擴(kuò)散、不同相間的物質(zhì)交換、相變、化學(xué)反應(yīng))取決于含氮物質(zhì)的各種熱力學(xué)與動(dòng)力學(xué)性質(zhì)等。這些過(guò)程也伴隨著能量的交換、轉(zhuǎn)化和力的傳遞。要想正確地理解有流體參與的地質(zhì)過(guò)程, 不但需要計(jì)算流體的性質(zhì), 還需要計(jì)算流體的膨脹系數(shù)()、壓縮系數(shù)()、聲速()、等壓或等容摩爾熱容(C,C)、焓()、熵()、Gibbs自由能()、Helmholtz自由能()、內(nèi)能()、化學(xué)勢(shì)()、逸度()、活度()等許多熱力學(xué)性質(zhì)。例如, 相平衡、化學(xué)平衡、熱效應(yīng)的計(jì)算則需要用到、、、、、等。另外, 一些流體動(dòng)力學(xué)性質(zhì)的計(jì)算也涉及熱力學(xué)性質(zhì)。例如, 流體的聲速和粘度()的計(jì)算需要用到流體的等壓及等容熱容、密度、膨脹系數(shù)和壓縮系數(shù)等。

上述的各種熱力學(xué)性質(zhì)都可以根據(jù)狀態(tài)方程、熱容、標(biāo)準(zhǔn)焓、標(biāo)準(zhǔn)熵以及相關(guān)的熱力學(xué)關(guān)系式推算出來(lái)。根據(jù)狀態(tài)方程, 很容易計(jì)算出膨脹系數(shù)、壓縮系數(shù)、逸度系數(shù)()、摩爾剩余焓(R)、摩爾剩余熵(R)等性質(zhì)。在此基礎(chǔ)上可以進(jìn)一步推算化學(xué)勢(shì)、逸度、活度等許多其他的熱力學(xué)性質(zhì)。實(shí)際上,、、R、R、或等性質(zhì)很少有實(shí)驗(yàn)或模擬數(shù)據(jù)的報(bào)道。不過(guò), 通過(guò)參考模型可以計(jì)算這些性質(zhì)。在此基礎(chǔ)上, 根據(jù)剩余性質(zhì)的定義及理想氣體的焓、熵與壓力的關(guān)系, 用參考模型算出真實(shí)流體的R和R, 然后很容易算出R和ln。

根據(jù)附錄A中的表達(dá)式和方程(17)預(yù)測(cè)ln、R、R和C–C。其中應(yīng)該采用方程(17)在給定溫壓下的體積根(附錄B)。這些性質(zhì)的預(yù)測(cè)結(jié)果與參考模型計(jì)算值的比較見(jiàn)圖12?？梢钥闯? 方程(17)與參考模型對(duì)ln和R的計(jì)算結(jié)果吻合得非常好。即使在低溫區(qū)(200～273.15 K), ln的計(jì)算結(jié)果也吻合得很好。R和C–C的計(jì)算結(jié)果大多數(shù)吻合得很好。即使C–C在275 K附近的幾個(gè)峰也被精確地重現(xiàn)。實(shí)際上, 由于狀態(tài)方程推導(dǎo)C–C的過(guò)程涉及與對(duì)溫度的二階導(dǎo)數(shù), 所以C–C對(duì)狀態(tài)方程的偏差很敏感, 其計(jì)算精度可以作為檢驗(yàn)狀態(tài)方程精度的一個(gè)很好的指標(biāo)。

(a)sim為MC模擬數(shù)據(jù), 據(jù)Etters et al. 1986; (b)sim為MD模擬數(shù)據(jù), 據(jù)Str?k and Krukowski, 2007; (c)sim為本研究的MD模擬數(shù)據(jù)。

圖12 方程(17)的預(yù)測(cè)結(jié)果與的參考模型計(jì)算值的比較(200～2000 K)(參考模型據(jù)Span et al., 2000)

不過(guò), 在低溫區(qū), 隨著壓力的增大,R和C–C的偏差逐漸增大。這主要是因?yàn)榉匠?17)的偏差隨著溫度的降低而逐漸增大。

此外,R在低溫、高壓區(qū)的偏差較大還有另外的因素。設(shè)Δln、ΔR與ΔR分別為ln、R與R計(jì)算結(jié)果的偏差。由于R=R–R,R=ln, 從而可得ΔR–ΔR=Δln。另一方面, 由于ln的計(jì)算結(jié)果很精確, 可以假設(shè)Δln≈0。由此可得ΔR≈ΔR/。因此, 當(dāng)ΔR相同時(shí), 溫度越低, ΔR越大。也就是說(shuō), 低溫放大了R的偏差。不過(guò), 這幾乎不影響Gibbs自由能計(jì)算結(jié)果的精確性。因?yàn)閘n的計(jì)算結(jié)果很精確, 所以R與R的偏差幾乎總能相互抵消。

7 結(jié) 論

為了精確地計(jì)算高溫、高壓下含N2地質(zhì)流體混合物的熱力學(xué)性質(zhì), 本研究根據(jù)N2的高精度數(shù)據(jù)建立了超臨界氮的一種高精度四參數(shù)立方型狀態(tài)方程。在273.15～4273.15 K、0～0.85 g·cm?3范圍內(nèi), 該方程的體積偏差在0.41%以?xún)?nèi), 平均體積偏差不到0.1%。新方程可以比較精確地應(yīng)用到5073.15 K或0.926 g·cm?3, 相應(yīng)的體積偏差在0.6%或1.0%以?xún)?nèi)。新方程可以較好地外延至200 K、1.03 g·cm?3, 相應(yīng)的體積偏差不超過(guò)2.5%。新方程也可近似地應(yīng)用到1.2 g·cm?3, 相應(yīng)的最高溫壓可達(dá)5573.15 K、7.5 GPa。

為了檢驗(yàn)新方程在高溫、高壓下的精度, 本文利用分子動(dòng)力學(xué)方法模擬了N2在1473.15～5073.15 K、0.3～6.0 GPa范圍內(nèi)的性質(zhì)。這些結(jié)果及現(xiàn)有的高溫或高壓實(shí)驗(yàn)數(shù)據(jù)、分子動(dòng)力學(xué)及Monte Carlo模擬數(shù)據(jù)(0.1～1.2 g·cm?3)均與新方程吻合得很好或較好, 相應(yīng)的偏差小于或接近于實(shí)驗(yàn)或模擬的不確定度。由新方程預(yù)測(cè)的膨脹系數(shù)、壓縮系數(shù)、逸度系數(shù)、剩余焓、剩余熵等熱力學(xué)性質(zhì)均與參考模型的計(jì)算值吻合得很好。這些結(jié)果表明, 新方程顯著地優(yōu)于現(xiàn)有的其他地質(zhì)流體狀態(tài)方程, 從而為高溫、高壓下多元含氮混合物熱力學(xué)性質(zhì)的精確模擬或預(yù)測(cè)提供了一個(gè)較好的純組分模型。

新方程簡(jiǎn)單、易用, 外延性也較好, 也易于向混合物推廣, 但在精確度上不可能超過(guò)參考模型。另外, 參考模型參數(shù)擬合所用的部分高溫或高壓實(shí)驗(yàn)數(shù)據(jù)存在一定的系統(tǒng)偏差或較大的不確定度, 其可靠性明顯不如中低溫壓下的實(shí)驗(yàn)數(shù)據(jù)。未來(lái)的高溫、高壓實(shí)驗(yàn)若能提供新的高精度數(shù)據(jù), 則新方程的參數(shù)還可以進(jìn)一步優(yōu)化。

附錄A 一些基本熱力學(xué)性質(zhì)的表達(dá)式

A.1 膨脹系數(shù)與壓縮系數(shù)

由于只是的隱函數(shù), (A1)式中第二個(gè)等式的變換采用了隱函數(shù)的求導(dǎo)法則?？紤]到文獻(xiàn)中一般都沒(méi)有與數(shù)據(jù), 本研究采用C–C來(lái)間接地檢驗(yàn)與計(jì)算結(jié)果的準(zhǔn)確性:

A.2 逸度系數(shù)

在計(jì)算流體組分的化學(xué)勢(shì)(即摩爾Gibbs自由能)之前, 必須先計(jì)算組分的逸度系數(shù):

在方程(17)中, 由于是和的顯函數(shù), 采用公式(A11)來(lái)推導(dǎo)ln更為方便:

式中=/(壓縮因子)。逸度系數(shù)及剩余焓、剩余熵的表達(dá)式需要對(duì)狀態(tài)方程進(jìn)行積分才能得到。積分結(jié)果的形式取決于方程(17)中二次項(xiàng)分母(+)的根的情況。(+)的根的判別式Δ=2。由于Δ可能因溫度的變化而改變符號(hào), 所以Δ的值有兩種情況: 當(dāng)≠0時(shí), Δ>0; 當(dāng)=0時(shí), Δ=0。當(dāng)Δ>0時(shí), 將方程(17)代入公式(A11)可得:

當(dāng)Δ=0時(shí),

A.3 摩爾剩余焓與摩爾剩余熵

對(duì)于顯壓型方程[=(,)]可以直接由狀態(tài)方程推導(dǎo)摩爾剩余內(nèi)能R(即摩爾剩余熱力學(xué)能)。然后再由公式(A15)算出摩爾剩余焓R:

若已知R, 可由公式(A15)～公式(A17)式算出R。因此, 下面只需考慮R的計(jì)算。

當(dāng)Δ>0時(shí), 將方程(17)代入公式(A16), 積分后可得

當(dāng)Δ=0時(shí), 將方程(17)代入公式(A16), 積分后即可得到R的表達(dá)式。當(dāng)Δ=0時(shí),=0。據(jù)此可將R的表達(dá)式簡(jiǎn)化為

當(dāng)Δ=0時(shí), 溫度只是一個(gè)特定的值(溫度軸上的一個(gè)點(diǎn))。因此可以說(shuō), 出現(xiàn)Δ=0的概率無(wú)窮小。在具體的計(jì)算中實(shí)際上無(wú)需考慮這種情況。

附錄B 方程(17)在給定-條件下的求解

首先考慮解析法。任一立方型狀態(tài)方程都可統(tǒng)一轉(zhuǎn)化成如下形式:

其中可以是壓縮因子()、摩爾體積()或數(shù)密度(即的倒數(shù))。然后, 將系數(shù)1～3代入求根公式即可得到體積根, 詳見(jiàn)郭濤等(2016)的附錄A。

其次, 也可以采用二分法或其他方法。其中二分法是求解各種非線性方程最常用的方法之一。對(duì)于一定溫度下的超臨界流體, 由狀態(tài)方程(17)算出的關(guān)系是一條單調(diào)的曲線。因此, 方程(17)在其適用范圍內(nèi)的任意-條件下只有唯一的體積根, 滿(mǎn)足二分法的應(yīng)用條件。

在給定-條件下, 應(yīng)用二分法求體積根的思路大體如下:

①首先定義一個(gè)誤差函數(shù)(,), 如(,)=cal(,)/1。其中cal是由方程(17)算出的壓力,為給定的壓力,是一個(gè)任意的摩爾體積。與此同時(shí), 為迭代計(jì)算設(shè)定一個(gè)允許的最大誤差, 如=1.0×10?8。

②為待求根設(shè)置兩個(gè)合適的初值(即1和2), 以便粗略地圈定待求根所在的區(qū)間[1,2]。如果取初值時(shí)沒(méi)有任何參考值可用, 可取1=,2=(1+)。其中為方程(17)的協(xié)體積,為一個(gè)極小的正值;可取0～1之間某個(gè)較小的正值, 如0.1或0.2。如果有實(shí)驗(yàn)值或其他值可作為參考(ref), 可令1=(1–)ref,2=(1–)ref, 其中ref為參考值。

③將溫度與初值1和2代入方程(17), 得到1=cal(,1),2=cal(,2), 然后算出兩個(gè)壓力下的偏差, 即1=(,1)＝1/1,2=(,2)=2/?1。

④若1·2>0, 說(shuō)明方程(17)在[1,2]內(nèi)無(wú)根, 然后返回步驟②, 調(diào)整初值, 使其滿(mǎn)足1·2<0; 若1·2<0, 說(shuō)明方程(17)在[1,2]內(nèi)有根, 然后進(jìn)入下一步的迭代循環(huán)。

⑤令3=(1+2)/2, 算出相應(yīng)的壓力3=cal(,3); 同時(shí)算出相應(yīng)的壓力偏差3=(,3)=3/?1。如果|3|<, 則3就是允許誤差范圍內(nèi)待求根的數(shù)值解。此時(shí)可直接跳到步驟⑧; 否則, 執(zhí)行下一步運(yùn)算。

⑥如果|3|>, 則更新[1,2]的邊界。具體方法: 若1·3>0, 即=1和=3時(shí)的壓力偏差同號(hào)。這說(shuō)明3比1更接近于待求根。此時(shí)應(yīng)將邊界值1更新, 即令1=3,1=3,1=3。若1·3<0, 則進(jìn)行相反的運(yùn)算, 即令1=2,1=2,1=2。

⑦重復(fù)步驟⑤和⑥。

⑧記錄迭代結(jié)果, 結(jié)束求根。

致謝:本文作者對(duì)清華大學(xué)化學(xué)工程系于養(yǎng)信教授在分子動(dòng)力學(xué)模擬方面的指導(dǎo)和大力支持深表感謝, 同時(shí)感謝兩位匿名審稿專(zhuān)家提出的建設(shè)性修改意見(jiàn)。

段振豪. 2010. 地質(zhì)流體狀態(tài)方程. 中國(guó)科學(xué): 地球科學(xué), 40(4): 393–413.

郭磊. 2017. 高溫高壓下二氧化碳狀態(tài)方程的評(píng)價(jià). 石家莊: 河北地質(zhì)大學(xué)碩士學(xué)位論文: 1–50.

郭濤, 胡家文, 王向輝, 靳麗花, 王艷哲. 2016. 超臨界水的熱力學(xué)模擬: 一種可用到4273 K、2 GPa的立方型狀態(tài)方程. 巖石學(xué)報(bào), 32(7): 2196–2208.

胡家文. 2002. 維里型和立方型地質(zhì)流體狀態(tài)方程的理論研究. 成都: 成都理工大學(xué)博士學(xué)位論文: 1–101.

黃榮榮. 1995. Burnett法--測(cè)定裝置. 江蘇石油化工學(xué)院學(xué)報(bào), 7(3): 8–13.

李旭平, 王晗, 孔凡梅. 2019. 高溫–超高溫變質(zhì)作用成因研究——來(lái)自華北克拉通西部孔茲巖帶和南非Kaapvaal克拉通西南部Namaqua活動(dòng)帶與Bushveld變質(zhì)雜巖體的啟示. 巖石學(xué)報(bào), 35(2): 295–311.

劉斌. 2011. 簡(jiǎn)單體系水溶液包裹體pH和Eh的計(jì)算. 巖石學(xué)報(bào), 27(5): 1533–1542.

劉守偈, 李江海, Santosh M. 2008. 內(nèi)蒙古土貴烏拉孔茲巖帶超高溫變質(zhì)作用: 變質(zhì)反應(yīng)結(jié)構(gòu)及指示. 巖石學(xué)報(bào), 24(6): 1185–1192.

盧煥章, 范宏瑞, 倪培, 歐光習(xí), 沈昆, 張文淮. 2004. 流體包裹體. 北京: 科學(xué)出版社: 1–485.

倪培, 遲哲, 潘君屹, 王國(guó)光, 陳輝, 丁俊英. 2018. 熱液礦床的成礦流體與成礦機(jī)制——以中國(guó)若干典型礦床為例. 礦物巖石地球化學(xué)通報(bào), 37(3): 370–394.

孫賀, 肖益林. 2009. 流體包裹體研究: 進(jìn)展、地質(zhì)應(yīng)用及展望. 地球科學(xué)進(jìn)展, 24(10): 1105–1120.

王艷哲, 2015. 高溫高壓下水的狀態(tài)方程的檢驗(yàn)與評(píng)價(jià). 石家莊: 石家莊經(jīng)濟(jì)學(xué)院碩士學(xué)位論文: 1–71.

魏春景. 2016. 麻粒巖相變質(zhì)作用與花崗巖成因–Ⅱ: 變質(zhì)泥質(zhì)巖高溫–超高溫變質(zhì)相平衡與S型花崗巖成因的定量模擬. 巖石學(xué)報(bào), 32(6): 1625–1643.

肖益林, 余成龍, 王洋洋, 陸一敢. 2018. 變質(zhì)作用與流體包裹體: 進(jìn)展與展望. 礦物巖石地球化學(xué)通報(bào), 37(3): 424–440.

徐航, 蔡旭東, 胡芃, 金熠. 2015. 基于磁懸浮密度計(jì)的高精度流體測(cè)量系統(tǒng). 中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào), 45(5): 404–408.

楊光璐, 藺玉秋, 劉加林, 鮑君剛. 2004. 遼河油田稠油油藏氮?dú)馀菽?qū)適應(yīng)性研究. 新疆石油地質(zhì), 22(2): 188–190.

翟偉, 孫曉明, 徐莉, 張澤明, 梁金龍, 梁業(yè)恒, 沈昆. 2005. 蘇北青龍山超高壓變質(zhì)榴輝巖流體包裹體特征與流體演化. 巖石學(xué)報(bào), 21(2): 482–488.

張澤明, 沈昆, 趙旭東, 石超. 2006. 超高壓變質(zhì)作用過(guò)程中的流體——來(lái)自蘇魯超高壓變質(zhì)巖巖石學(xué)、氧同位素和流體包裹體研究的限定. 巖石學(xué)報(bào), 22(7): 1985–1998.

周健, 陸小華, 王延儒, 時(shí)鈞. 1999. 超臨界水的分子動(dòng)力學(xué)模擬. 物理化學(xué)學(xué)報(bào), 15(11): 1017–1122.

Andersen T, Austrheim H, Burke E A J, Elvevold S. 1993. N2and CO2in deep crustal fluids: Evidence from the Caledonides of Norway., 108(1–4): 113–132.

Andersen T, Burke E A J, Neumann E-R. 1995. Nitrogen-rich fluid in the upper mantle: Fluid inclusions in spinel dunite from Lanzarote, Canary Islands., 120(1): 20–28.

Antanovich A A, Plotnikov M A. 1976. Experimental determination of the density of nitrogen at high pressures and temperatures., 21(2): 99–100.

Antanovich A A, Plotnikov M A. 1977. Investigation of the thermodynamic properties of nitrogen at pressures up to 8 kbar and temperatures up to 1800 °K., 33(2): 929–934.

Artimenko M V. 2017. Two-phase fluid induced by N2in metamorphic rocks., 475: 40–51.

Atilhan M. 2007. High accuracymeasurements up to 200 MPa between 200 K and 500 K using a compact single sinker magnetic suspension densimeter for pure and natural gas like mixtures. College Station: Texas Agricultural and Mechanical University: 1–141.

Bakker R J. 1999. Adaptation of the Bowers and Helgeson (1983) equation of state to the H2O-CO2-CH4-N2-NaCl system., 154(1–4): 225–236.

Bakker R J. 2003. Package1. Computer programs for analysis of fluid inclusion data and for modelling bulk fluid properties., 194(1–3): 3–23.

Belak J, Etters R D, Lesar R A. 1988. Thermodynamic properties and equation of state of dense fluid nitrogen.89(3): 1625–1633.

Belonoshko A B, Saxena S K. 1992. A unified equation of state for fluids of C-H-O-N-S-Ar composition and their mixtures up to very high temperatures and pressures., 56(10): 3611–3626.

Berkesi M, Káldos R, Park M, Szabó C, Váczi T, T?r?k K, Nemeth B, Czuppon G. 2017. Detection of small amounts of N2in CO2-rich high-density fluid inclusions from mantle xenoliths., 29(3): 423–431.

Boates B, Bonev S A. 2009. First-order liquid-liquid phase transition in compressed nitrogen., 102(1), 015701.

Brennan J K, Rice B M. 2002. Molecular simulation of shocked materials using the reactive Monte Carlo method., 66(2), 021105.

Chamorro C R, Segovia J J, Martín M C, Villama?án M A, Estela-Uribe J F, Trusler J P M. 2006. Measurement of the (pressure, density, temperature) relation of two (methane + nitrogen) gas mixtures at temperatures between 240 and 400 K and pressures up to 20 MPa using an accurate single-sinker densimeter., 38(7): 916–922.

Chapman D S. 1986. Thermal gradients in the continental crust.,,, 24(1): 63–70.

Chen Q, Zhang Z, Wang Z, Li W C, Gao X Y, Ni H. 2019.Raman spectroscopic study of nitrogen speciation in aqueous fluids under pressure., 506: 51–57.

Cheng S Y, Shang F, Ma W G, Jin H, Sakoda N, Zhang X, Guo L J. 2019. Density measurements of the H2-CO2-CH4-CO-H2O system by the isochoric method at 722?930 K and 15.4?30.3 MPa., 64(9): 4024–4036.

Christoforakos M, Franck E U. 1986. An equation of state for binary fluid mixtures to high temperatures and high pressures., 90(9): 780–789.

Churakov S V, Gottschalk M. 2003. Perturbation theory based equation of state for polar molecular fluids: I. Pure fluids., 67(13): 2397–2414.

Dipietro J A, 2013. Landscape Evolution in the United States. New York: Elsevier: 327–344.

Driver K P, Militzer B. 2016. First-principles equation of state calculations of warm dense nitrogen.93(6), 064101.

Duan Z H, MФller N, Weare J H. 1992. Molecular dynamics simulation ofproperties of geological fluids and a general equation of state of nonpolar and weakly polar gases up to 2000 K and 20, 000 bar., 56(10): 3839–3845.

Duan Z H, MФller N, Weare J H. 1996. A general equation of state for supercritical fluids and molecular dynamics simulation of mixtureproperties., 60(7): 1209–1216.

Dubessy J, Poty B, Ramboz C. 1989. Advances in C-O-H-N-Sfluid geochemistry based on micro-Raman spectrometricanalysis of fluid inclusions., 1(4): 517–534.

Dubois V, Desbiens N, Clérouin J. 2017. Seeking an accurate generalized-gradient approximation functional for high pressure molecular fluids., 122(18), 185902.

Elkins L J, Fischer T P, Hilton D R, Sharp Z D, McKnight S, Walker J. 2006. Tracing nitrogen in volcanic and geothermalvolatiles from the Nicaraguan volcanic front., 70(20): 5215–5235.

Elvevold S, Andersen T. 1993. Fluid evolution during metamorphism at increasing pressure: Carbonic- and nitrogen-bearing fluid inclusions in granulites from Oksfjord, north Norwegian Caledonides., 114(2): 236–246.

Etters R D, Belak J, LeSar R. 1986. Thermodynamic character of the vibron frequencies and equation of state in dense, high-temperature, fluid N2., 34(6): 4221–4223.

Evers C, L?sch H W, Wagner W. 2002. An absolute viscometer-densimeter and measurements of the viscosity of nitrogen, methane, helium, neon, argon, and krypton over a wide range of density and temperature., 23(6): 1411–1439.

Fischer T P. 2008. Fluxes of volatiles (H2O, CO2, N2, Cl, F) from arc volcanoes.42(1): 21–38.

Fu B, Touret J L R, Zheng Y F. 2001. Fluid inclusions in coesite-bearing eclogites and jadeite quartzite at Shuanghe, Dabie Shan (China)., 19(5): 531–547.

Fu Z J, Chen Q F, Li Z G, Tang J, Zhang W, Quan W L, Li J T, Zheng J, Gu Y J. 2019.study of the structures and transport properties of warm dense nitrogen., 31: 52–58.

García-Jarana M B, Kings I, Sánchez-Oneto J, Portela J R, Al-Duri B. 2013. Supercritical water oxidation of nitrogen compoundswith multi-injection of oxygen., 80: 23–29.

Gottschalk M. 2007. Equations of state for complex fluids.65(1): 49–97.

Heilig M, Franck E U. 1989. Calculation of thermodynamic properties of binary fluid mixtures to high temperatures and high pressures., 93(8): 898–905.

Holland T, Powell R. 1991. A Compensated-Redlich-Kwong (CORK) equation for volumes and fugacities of CO2and H2O in the range 1 bar to 50 kbar and 100–1600 ℃., 109(2): 265–273.

Holloway J R, 1977. Fugacity and activity of molecular species in supercritical fluids // Fraser D G (Ed.). Thermodynamics in Geology. Dordrecht: D. Reidel Publishing Company: 161–181.

Johnson R T, Stewart R B, Jahangiri M. 1986. Thermodynamic properties of nitrogen from the freezing line to 2000 K at pressures to 1000 MPa., 15(2): 735–908.

Johnson B, Goldblatt C. 2015. The nitrogen budget of Earth.148: 150–173.

Johnson J D, Shaw M S, Holian B L. 1984. The thermodynamics of dense fluid nitrogen by molecular dynamics., 80(3): 1279–1294.

Kraussler M, Binder M, Fail S, Bosch K, Hackel M, Hofbauer H. 2016. Performance of a water gas shift pilot plant processing product gas from an industrial scale biomass steam gasification plant., 89: 50–57.

Kress J D, Mazevet S, Collins L A, Wood W W. 2000. Density-functional calculation of the Hugoniot of shocked liquid nitrogen., 63(2), 024203.

Krooss B M, Littke R, Müller B, Frielingsdorf J, Schwochau K, Idiz E F. 1995. Generation of nitrogen and methane from sedimentary organic matter: Implications on the dynamics of natural gas accumulations., 126(3–4): 291–318.

Krukowski S, Str?k P. 2006. Equation of state of nitrogen (N2) at high pressures and high temperatures: Moleculardynamics simulation., 124(13), 134501.

Lambert R S J, 1983. Metamorphism and thermal gradients in the Proterozoic continental crust // Medaris L G J, Byers C W, Mickelson D M, Shanks W C (Eds.). ProterozoicGeology: Selected Papers from an International Proterozoic Symposium. Boulder: Geological Society of America: 155–165.

Lee B I, Kesler M G. 1975. A generalized thermodynamic correlation based on three-parameter corresponding states., 21(3): 510–527.

Lv L, Zhang L, Yang M L. 2018. Understanding the phase separation of N2/H2O and CO2/H2O binary systems through reactive force fields-based molecular dynamics simulations., 124(23), 235901.

Maiti A, Gee R H, Bastea S, Fried L E. 2007. Phase separation in H2O-N2mixture: Molecular dynamics simulations using atomistic force fields., 126(4), 044510.

Malbrunot P, Vodar B. 1973. Experimentaldata and thermodynamic properties of nitrogen up to 1000 ℃ and 5000 bar., 66(2): 351–363.

Mallik A, Li Y, Wiedenbeck M. 2018. Nitrogen evolution within the Earth’s atmosphere-mantle system assessed by recycling in subduction zones., 482: 556–566.

Mantilla I D, Cristancho D E, Ejaz S, Hall K R. 2010. Newdata for nitrogen at temperatures from (265 to 400) K at pressures up to 150 MPa., 55(10): 4227–4230.

Mao S D, Lu M X, Shi Z M. 2017. Prediction of theand VLE properties of natural gases with a general Helmholtz equation of state. Part I: Application to the CH4-C2H6-C3H8-CO2-N2system., 219: 74–95.

Marko F, Hurai V, Dyda M, Almeida G, Prochaska W, Thomas R. 2006. Tectonic and fluid inclusion constraints on the origin of quartz veins with giant crystals in the Tocantins structural province (Cristalandia, central Brazil)., 21(3): 239–251.

Mazevet S, Johnson J D, Kress J D, Collins L A, Blottiau P. 2001. Density-functional calculation of multiple-shock Hugoniots of liquid nitrogen., 65(1), 014204.

Mikhail S, Barry P H, Sverjensky D A. 2017. The relationship between mantle pH and the deep nitrogen cycle., 209: 149–160.

Mills R L, Liebenberg D H, Bronson J C. 1975. Sound velocity and the equation of state of N2to 22 kbar., 63(3): 1198–1204.

Mochalov M A, Zhernokletov M V, Il’kaev R I, Mikhailov A L, Fortov V E, Gryaznov V K, Iosilevskiy I L, Mezhevov A B, Kovalev A E, Kirshanov S I, Grigor’eva Y A, Novikov M G, Shuikin A N. 2010. Measurement of density, temperature, and electrical conductivity of a shock-compressed nonideal nitrogen plasma in the megabar pressure range., 110(1): 67–80.

Mourato M F B, Calado J C G, Palavra A M F. 1999. Automated isochoric apparatus for pressure, density, temperature measurements of binary gaseous mixsures at high temperatures., 31(1): 91–98.

Nellis W J, Radousky H B, Hamilton D C, Mitchell A C, Holmes N C, Christianson K B, van Thiel M. 1991. Equation-of-state, shock-temperature, and electrical- conductivity data of dense fluid nitrogen in the region of the dissociative phase transition., 94(3): 2244–2257.

Patil P V. 2005. Commissioning of a magentic suspension densitometer for high-accuracy density measurements of natural gas mixtures. College Station: Texas A&M University: 1–331.

Perkins R A, Roder H M, Decastro C A N. 1991. A high-temperature transient hot-wire thermal conductivityapparatus for fluids., 96(3): 247–269.

Radousky H B, Ross M. 1988. Shock-induced cooling in high-density fluid nitrogen., 1(1): 39–52.

Redlich O, Kwong J N S. 1949. On the thermodynamics of solutions, V. An equation of state. Fugacities of gaseous solutions., 44(1): 233–244.

Robertson S L, Babb S E. 1969. Isotherms of nitrogen to 400 ℃and 10000 bar., 50(10): 4560–4564.

Roder H M. 1981. A transient hot wire thermal conductivity apparatus for fluids., 86(5): 457–493.

Roskosz M, Bouhifd M A, Jephcoat A P, Marty B, Mysen B O. 2013. Nitrogen solubility in molten metal and silicate at high pressure and temperature., 121: 15–28.

Sakoda N, Shindo K, Motomura K, Shinzato K, Kohno M, Takata Y, Fujii M. 2012. Burnett method with absolute pressure transducer and measurements forproperties of nitrogen and hydrogen up to 473 K and 100 MPa., 33(1): 6–21.

Saleh B, Wendland M. 2005. Measurement of vapor pressures and saturated liquid densities of pure fluids with a new apparatus., 50(2): 429–437.

Saxena S, Fei Y. 1987a. Fluids at crustal pressures and temperatures I. Pure species., 95(3): 370–375.

Saxena S, Fei Y. 1987b. High pressure and high temperature fluid fugacities., 51(4): 783–791.

Siemann M G, Ellendorff B. 2001. The composition of gases in fluid inclusions of late Permian (Zechstein) marine evaporites in Northern Germany., 173(1–3): 31–44.

Smith E M, Kopylova M G, Frezzotti M L, Afanasiev V P. 2014. N-rich fluid inclusions in octahedrally-grown diamond., 393: 39–48.

Soave G S. 1995. A noncubic equation of state for the treatment of hydrocarbon fluids at reservoir conditions., 34(11): 3981–3994.

Soave G S. 1999. An effective modification of the Benedict- Webb-Rubin equation of state., 164(2): 157–172.

Sobolev N V, Logvinova A M, Tomilenko A A, Wirth R, Bul’bak T A, Luk'yanova L I, Fedorova E N, Reutsky V N, Efimova E S. 2019. Mineral and fluid inclusions in diamonds from the Urals placers, Russia: Evidence for solid molecular N2and hydrocarbons in fluid inclusions., 266: 197–219.

Sokol A G, Palyanov Y N, Tomilenko A A, Bul’bak T A, Palyanova G A. 2017. Carbon and nitrogen speciation in nitrogen-rich C–O–H–N fluids at 5.5–7.8 GPa., 460: 234–243.

So?nicka M, Lüders V. 2019. Fluid inclusion evidence for low-temperature thermochemical sulfate reduction (TSR) of dry coal gas in Upper Permian carbonate reservoirs (Zechstein, Ca2) in the North German Basin., 534: 119453.

Span R, Lemmon E W, Jacobsen R T, Wagner W, Yokozeki A. 2000. A reference equation of state for thethermodynamic properties of nitrogen for temperatures from 63.151 to 1000 K and pressures to 2200 MPa., 29(6): 1361–1433.

Str?k P, Krukowski S. 2007. Molecular nitrogen-N2properties: The intermolecular potential and the equation of state., 126(19): 194501.

Sun L, Ely J F. 2004. Universal equation of state for engineering application: Algorithm and application to non-polar and polar fluids., 222–223: 107–118.

Sun R, Lai S, Dubessy J. 2015. Calculations of vapor-liquid equilibria of the H2O-N2and H2O-H2systems with improved SAFT-LJ EOS., 390: 23–33.

Tarney J, Windley B F. 1977. Chemistry, thermal gradients and evolution of the lower continental crust., 134(2): 153–172.

Tsiklis D S, Polyakov E V. 1968. Measuring the compressibility of gases by the displacement method. nitrogen compressibilityat pressures up to 10,000 atm and temperatures to 400°., 12(9): 901–904.

Van den Kerkhof A M, Touret J L R, Maijer C, Jansen J B H. 1991. Retrograde methane-dominated fluid inclusions from high-temperature granulites of Rogaland, southwestern Norway., 55(9): 2533–2544.

Wang Y C, Wang K Y, Konare Y. 2018. N2-rich fluid in the vein-type Yangjingou scheelite deposit, Yanbian, NE China., 8(1), 5662.

Xiao Y, Hoefs J, van den Kerkhof A M, Simon K, Fiebig J, Zheng Y F. 2002. Fluid evoluition during HP and UHP metamorphism in Dabie Shan, China: Constaints from mineral chemistry, fluid inclusions and stable isotopes., 43(8): 1505–1527.

Yang J, Ren Y, Tian A M, Sun H. 2000. COMPASS force field for 14 inorganic molecules, He, Ne, Ar, Kr, Xe, H2, O2, N2, NO, CO, CO2, NO2, CS2, and SO2, in liquid phases. Journal of Physical Chemistry B, 104(20) : 4951–4957.

Yang X Q, Li Z L, Yu S Q. 2016. Phase equilibrium modeling, fluid inclusions and origin of charnockites in the Datian region of the northeastern Cathaysia Block, South China., 126: 14–28.

Yang X X, Richter M, Wang Z, Li Z. 2015. Density measurements on binary mixtures (nitrogen + carbon dioxide and argon + carbon dioxide) at temperatures from (298.15 to 423.15) K with pressures from (11 to 31) MPa using a single-sinker densimeter.Chemical Thermodynamics, 91: 17–29.

Ye L, Liu Y, Yang Y, Gao W. 2012. The characteristics and significance of pure nitrogen fluid inclusions in Xikuangshan copper deposit, Dongchuan, Yunnan of China., 31(1): 78–84.

Yoshioka T, Wiedenbeck M, Shcheka S, Keppler H. 2018. Nitrogen solubility in the deep mantle and the origin of Earth’s primordial nitrogen budget., 488: 134–143.

Zhang Z G, Zhang C, Geng M. 2016. Equations of state for aqueous solutions under mantle conditions., 59(6): 1095–1106.

Four-parameter cubic equation of state and molecular dynamics simulation of supercritical nitrogen

JIANG Siyu1, 2, GUO Tao1, 3, HU Jiawen1, 2*

(1. Hebei Key Laboratory of Strategic Critical Mineral Resources, Hebei GEO University, Shijiazhuang 050031, Hebei, China; 2. College of Earth Science, Hebei GEO University, Shijiazhuang 050031, Hebei, China; 3. School of Mathematics and Science, Hebei GEO University, Shijiazhuang 050031, Hebei, China)

A four-parameter cubic equation of state for supercritical nitrogen (N2) was developed using highly accurate pressure-volume-temperature () data and their extrapolation results. The conditions that best fit the new equation are 273–4273 K and 0–0.85 g·cm–3. In this range, the average and maximum volume deviations of the new equation are less than 0.1% and 0.41%, respectively; the range of pressure of an isotherm increases with increasing temperature, so the maximum pressure may rise to 2.9 GPa. When density increases to 1.03 g·cm–3and temperature increases to 5000 K, the volume deviations are within 0.6%–2.5%, and the maximum pressure increases to 3.3–5.0 GPa. The new equation can also be approximately used up to 1.2 g·cm–3, below which the maximum temperature and pressure reach 5573 K and 7.5 GPa, respectively. In addition, the PVT properties of N2at 1473–5073 K and 0.3–6.0 GPa were simulated using molecular dynamics. These results, along with the existing experimental data and molecular dynamics or Monte Carlo simulation results at high temperatures or densities (0.1–1.2 g·cm–3), are in good or excellent agreement with the new equation. The fugacity coefficient, residual enthalpies, residual entropies, and other thermodynamic properties predicted by the new equation agree very well with those calculated by the reference model. These results indicate that the new equation is significantly better than existing cubic, virial-type, and more complex equations of state for geofluids.

supercritical nitrogen; equation of state;properties; fugacity coefficient; residual enthalpy; molecular dynamics

O41; O642

A

0379-1726(2022)03-0344-21

10.19700/j.0379-1726.2022.03.008

2020-07-07;

2020-11-04

河北省自然科學(xué)基金(D2018403089)、河北省高等學(xué)校科學(xué)技術(shù)研究重點(diǎn)項(xiàng)目(ZD2018026)和河北省研究生創(chuàng)新資助項(xiàng)目(CXZZSS2020112)聯(lián)合資助。

姜思雨(1996–), 男, 碩士研究生, 地球化學(xué)專(zhuān)業(yè)。E-mail: jiang_geol@sina.com

胡家文(1967–), 男, 教授, 主要從事計(jì)算地球化學(xué)研究。E-mail: hu_jiawen@sina.com