超临界氮的四参数立方型状态方程与分子动力学模拟

姜思雨, 郭 涛, 胡家文*

超临界氮的四参数立方型状态方程与分子动力学模拟

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

(1. 河北地质大学 河北省战略性关键矿产资源研究重点实验室, 河北 石家庄 050031; 2. 河北地质大学 地球科学学院, 河北 石家庄 050031; 3. 河北地质大学数理教学部, 河北 石家庄 050031)

根据精确的压力–摩尔体积–温度()数据及其外延结果, 本研究建立了超临界氮(N2)的四参数立方型状态方程。新方程最适宜的条件为273～4273 K、0～0.85 g·cm–3。在此范围内, 新方程的最大体积偏差在0.41%以内, 平均体积偏差在0.1%以内; 等温线的压力范围随着温度的上升而扩大,最大压力达2.9 GPa。当密度增至1.03 g·cm–3、温度增至5000 K时, 体积偏差在0.6%～2.50%以内, 最高压力达3.3～5.0 GPa。新方程也可近似地应用到密度为1.2 g·cm–3。此时相应的温压可达5573 K、7.5 GPa。此外, 本研究还用分子动力学方法模拟了超临界氮在1473～5073 K、0.3～6.0 GPa范围内的性质。所得结果以及现有的高温或高密度实验数据、分子动力学和蒙特卡洛(Monte Carlo)模拟数据(0.1～1.2 g·cm–3)均与新方程吻合得很好或较好。由新方程预测的逸度系数、剩余焓、剩余熵及其他热力学性质与参考模型的计算值均很吻合。这些结果表明新方程显著地优于现有的立方型、维里型和更复杂的地质流体状态方程。

超临界氮; 状态方程;性质; 逸度系数; 剩余焓; 分子动力学

0 引 言

N2不仅是太阳系内巨行星及太阳系外行星的主要成分之一(Mochalov et al., 2010; Driver and Militzer, 2016; Dubois et al., 2017), 也是地球表面及内部常见的流体组分之一。由于N2的临界温度(C)很低(126.192 K)、临界压力(C)也不太高(3.3958 MPa), 因此N2在自然界经常以超临界流体形式存在。大量研究结果表明, 含N2流体可以出现在不同的地质环境中。例如, 含N2的流体包裹体常见于各种变质岩(Fu et al., 2001; 张泽明等, 2006; Wang et al., 2018; 肖益林等, 2018)、幔源金刚石(Smith et al., 2014; Sobolev et al., 2019)、海相蒸发岩(Siemann and Ellendorff, 2001)、地热流体(Dubessy et al., 1989)、煤层(Sośnicka and Lüders, 2019)和多种矿床(倪培等, 2018), 其中有的高密度包裹体在室温下甚至呈现固态(Sobolev et al., 2019)。这些包裹体保存了地质环境的许多信息, 通过对流体包裹体的组成、相变温度(冰点、水合物的形成或熔化温度、包裹体的部分或整体均一化温度)等性质的分析和相关的热力学模拟可以获取地质过程的温度、压力和流体密度等重要物理化学条件及其变化值(卢焕章等, 2004), 进而计算流体系统的氧逸度(O2)、pH值等(刘斌, 2011; Artimenko, 2017)。这些结果对于理解流体–矿物的相互作用、矿物、岩石和矿床(包括油气藏)的形成与演化、变质作用、地热与岩浆活动等许多过程十分重要(卢焕章等, 2004; 孙贺和肖益林, 2009; 倪培等, 2018; 肖益林等, 2018)。此外, 含N2的流体也会出现在煤层气(Krooss et al., 1995; Sośnicka and Lüders, 2019)、地幔捕虏体(Andersen et al., 1995; Berkesi et al., 2017)、熔体或岩浆(Yoshioka et al., 2018)、火山气体与地热流体(Elkins et al., 2006; Fischer, 2008)。地质流体中的N2可以以微量、少量或主量组分甚至纯N2的形式出现(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)。这些流体及其地质环境的温压条件可以从常温、常压变化到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)。

在地球或天体系统内, 分子态N2的主要存在形式是流体。此外, 氮还有许多其他存在形式, 包括无机、有机的离子或分子形式(固态或流体态)。其中, 分子态的流体氮与氮的其他形式之间通过迁移和反应相互转化(Roskosz et al., 2013; Johnson and Goldblatt, 2015; Sokol et al., 2017; Mallik et al., 2018)。因此, 要正确地理解含氮物种之间的转化(包括相平衡和化学平衡)及深部氮循环, 必须掌握含氮物质的热力学性质(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也是天然气的主要成分之一。天然气工业的许多环节都需要了解N2及相关混合物的热力学性质。N2是许多化工生产、超临界水氧化、废水处理、生物质加工、油田后期开采过程的原料或产物(周健等, 1999; 杨光璐等, 2004; García-Jarana et al., 2013; Kraussler et al., 2016), 也是许多含氮物质燃烧、爆炸或冲击压缩过程的主要产物(Driver and Militzer, 2016; Dubois et al., 2017; Fu et al., 2019)。

对上述各种过程的理解、模拟、设计和控制需要准确地掌握N2及其混合物的压力–摩尔体积–温度–组成()、热容、焓、熵及许多其他的热力学性质, 其中和热容是最基本的性质。和热容等基础数据被总结成各种热力学模型, 如状态方程、Helmholtz和Gibbs自由能模型。这些模型与微积分运算相结合可以推算出其他热力学性质。

1 N2的PVT实验或模拟数据

迄今为止, 人们已通过实验手段获得流体N2的大量热力学数据, 包括几组极高温压下冲击压缩实验的数据(Radousky and Ross, 1988; Nellis et al., 1991; Mochalov et al., 2010), 而其余实验数据的最高温压只有1800 K和2.2 GPa。Span et al. (2000)对截止到1997年的N2的各种热力学实验数据做了比较全面的搜集、统计和不确定度评估, 尽管也有一些实验数据被他们遗漏(Malbrunot and Vodar, 1973; Antanovich and Plotnikov, 1977; Roder, 1981; Belak et al., 1988; Perkins et al., 1991; 黄荣荣, 1995)。在近20年中, 陆续有一些N2的实验数据发表, 但均局限在中低温压范围(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)。总体而言, 绝大多数流体N2的热力学实验数据都在1000 K、1 GPa以下。其中相当一部分数据存在较大的不确定度或系统偏差; 而且,随着温度或压力的升高, 不确定度或系统偏差增加较快。对于许多高温高压生产或相关研究(如燃烧、含能材料、爆轰与冲击波物理、地球与行星深部的物理化学研究)而言, 这种情况远不能满足实际需要。

为了获取更多的高温高压热力学数据, 人们采用原子–分子尺度上的计算机模拟, 如分子动力学(molecular dynamics, MD)、蒙特卡洛(monte carlo, MC)、密度泛函理论(density function theory, DFT)、第一性原理(first principles, FP)模拟、量子力学从头计算法(), 或这些方法的某种联合应用。这些方法只需少量参数或者无需任何经验参数就可以获得很宽温压范围内(包括极端温压下)流体的单相、多相系统和化学反应的各种热力学数据, 从而极大地拓宽了流体热力学数据的温压范围, 同时也显示了很好的预测能力(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)。与传统的高温、高压实验方法相比, 这类方法在时间、成本、安全性和适用范围等方面具有巨大的优势。

目前, 通过实验(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的热力学数据。对于地球与行星深部研究而言, 1000～5000 K、1～5 GPa范围内流体N2的热力学数据特别重要, 但相应的实验数据却很少。相关的模拟数据虽然较多一些, 但有不少数据具有显著的系统偏差(Johnson et al., 1984; Driver and Militzer, 2016), 而且很多模拟数据只有图形报道(Etters et al., 1986; Strąk and Krukowski, 2007), 难以精确地再现(Maiti et al., 2007)或报道有误(Krukowski and Strąk, 2006)。因此, 目前仍然需要补充高温、高压下N2的热力学数据。

2 N2的热力学模型

自20世纪70年代以来, 许多研究致力于发展地质流体(包括含N2流体)的热力学模型, 其中主要是状态方程(Bakker, 2003; Gottschalk, 2007; 段振豪, 2010; Zhang et al., 2016)。目前, 含N2流体的热力学模型大致可以分为4类:

①立方型状态方程。自从van der Waals (vdW)提出第一个立方型状态方程以来, 不断有新的立方型状态方程发表。Redlich and Kwong (1949)(RK49)将vdW方程的引力参数改为/1/2, 所得方程如下:

式中:为协体积(co-volume);为通用气体常数;为摩尔体积。此方程在中低压(几十MPa)下的精度得到了显著的提高, 但在中高压下仍然欠佳。尽管如此, Holloway (1977)直接将该方程用于除H2O和CO2以外的非极性或弱极性地质流体组分(包括N2)。为了计算含CH4和N2的地质流体的热力学性质, Bakker (1999)(B99)将RK49方程的参数改为温度的线性函数:

式中:0和1为常量;0= 273.15 K。此方程中和单位分别为K、MPa和cm3·mol–1。对于N2,= 27 cm3·mol−1。Bakker将此方程应用到1300 K、1 GPa, 但偏差较大。

总体而言, 这类方程形式简单, 参数较少, 在给定温压之下有解析形式的体积根, 一般都具有较好的预测性, 便于以同样的形式推广到多组分混合物, 而且所需要的混合参数较少, 应用方便。目前, 这类方程数以百计, 在各种工业流体的热力学计算或模拟中得到了广泛的应用。不过, 立方型方程所适用的温压范围一般都在几百K和100 MPa以内, 在中高密度条件下对摩尔体积及相关广延性质的计算偏差通常较大或很大。

②维里型状态方程及其衍生形式。此类方程也较多, 但其中很少有方程能在很宽的–范围内保持较高的精度。例如, Lee and Kesler (1975)和Soave (1995, 1999)方程主要面向工业流体, 不适合高压地质流体。Saxena and Fei (1987a, 1987b)提出两种对应态形式的维里型方程, 其中前者Saxena and Fei (1987a)采用了分段函数, 后者Saxena and Fei (1987b)在中低压下不正确。但它们的偏差都较大或很大。Holland and Powell (1991)(HP91)将RK49方程近似地转化为以和为自变量的分式方程, 并加上2个补偿项:

式中:、和的单位分别为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。不过, 他们并未报道其方程对N2的具体效果。

Duan et al. (1992)(DMW92)利用CH4的实验数据和MD模拟结果, 针对非极性的地质流体组分(包括N2)建立了一种对应态形式的通用维里型状态方程:

Belonoshko and Saxena (1992)(BS92)利用MD模拟结果和exp-6势能参数建立了一种对应态形式的维里型状态方程:

③半经验的统计力学状态方程。此类方程中有一些不太复杂的方程(Christoforakos and Franck, 1986; Heilig and Franck, 1989), 但多数都比较复杂或很复杂, 如目前被广泛应用的有统计缔合流体理论(statistical associating fluid theory, SAFT)、微扰理论(perturbation theory, PT)状态方程(Churakov and Gottschalk, 2003; Sun et al., 2015), 其中SAFT方程一般所需参数较少, 预测性较好, 但对密度或体积的计算偏差较大。这类方程主要用于工业流体, 相应的温压范围比较有限, 一般不适于高温、高压下的地质流体。

④Helmholtz自由能模型。其中有一类是比较复杂的专用模型, 如Span et al. (2000)和Jacobsen et al. (1986)的Helmholtz自由能模型。这类模型适用温压范围较宽, 也很精确, 但需要拟合的参数很多, 很难以同样的形式推广到混合物。另一类是经过简化、具有一定通用性的Helmholtz自由能模型, 例如Sun and Ely (2004)与Mao et al. (2017)的模型。这类模型比较精确, 但比维里型方程要稍复杂一些, 能够应用的温压范围比较有限, 不适合于高温、高压地质流体。

关于含N2地质流体的上述热力学模型一般均有较大的偏差, 而且实际可用的温压范围与原作者所声明的结果经常有很大(或较大)的差别。这也印证了本研究组近期关于地质流体状态方程的研究结果(王艳哲, 2015; 郭涛等, 2016; 郭磊, 2017)。

目前, 关于含N2地质流体的热力学模型虽然很多, 但简单、精确、适用范围宽广的热力学模型仍然比较缺乏。鉴于目前地质流体组分及混合物的高温、高压热力学数据仍然短缺, 采用比较简单的模型来模拟地质流体的热力学性质是一种明智的选择。

3 四参数立方型状态方程的建立

考虑到立方型状态方程具有多方面的优点, 在适用范围的扩展方面仍有很大的潜力, 本研究拟建立一种新的立方型状态方程, 力图在较宽的温压范围内精确地再现或预测含N2流体的多种热力学性质。为此, 首先需要建立纯N2流体的状态方程。

3.1 方程的表达式

为了兼顾简洁性、精确度和适用范围, 本文采用以下的四参数立方型状态方程:

式中:2、3和均为温度的函数;为常数;= 8.31451 J·K–1·mol–1;、与的单位分别是MPa、cm3·mol−1和K。在Lennard-Jones流体第二维里系数的统计力学解析表达式中, 温度的指数均为1/4的整数倍。据此, 胡家文(2002)采用温度指数为1/4的整数倍(或相近值)的简单多项式来拟合第二维里系数, 效果很好。郭涛等(2016)关于水的立方型状态方程中引力参数采用了类似的形式。结果表明, 这种表达式可以有效地外延到极高的温度。

对方程(17)中2、3与的关系, 本研究尝试了不同的表达式(包括改变表达式的形式或增减参数的数量), 结果发现以下的形式效果最好:

3.2 方程的参数拟合

Span et al. (2000)从截止到1997年的N2的各种热力学数据中选用最好的数据建立了一种高精度的Helmholtz自由能模型。该模型最适宜的温压应用范围是63.151～1000 K、0～2.2 GPa。从三相点到523 K、0～12 MPa或240～523 K、0～30 MPa, 密度的不确定度为0.02%; 当压力很高(>1 GPa)时, 密度的不确定度可能会达到0.6%; 在270～350 K、12 MPa以下, 密度的不确定度仅为0.01%。鉴于该模型的精度和可靠性, 美国国家标准与技术研究所的官方网站将它推荐为N2的参考模型。不过, 该模型的形式很复杂, 线性和非线性参数都很多, 不便于拟合或应用, 也难以同样的形式应用于混合物, 因而其主要作用在于为科学研究提供可靠的数据, 有时也被用作混合物热力学模型中的纯组分模型。

为了保证新方程的准确性及外延结果的可靠性, 方程(17)参数拟合时所用的数据全部来自上述参考模型。为了使方程(17)具有尽可能宽的适用范围, 本研究对参考模型在上述范围内的计算结果向高温、高压方向进行了外延。

对于大多数流体而言, 中低密度时的等容线非常接近于直线; 在中高密度条件下, 等容线在中低温段会有明显的弯曲, 但在高温段弯曲程度会明显减弱。据此很容易对等容线进行精确的拟合, 同时也适于做大幅度的外延。本研究对N2等容线的拟合与外延均采用如下的多项式:

外延结果与参考模型在2000 K、2.2 GPa以下的计算结果一起用于方程(17)的参数拟合。

由于流体在临界点附近具有奇异性, 临界点附近性质的拟合难度远高于远离临界点的超临界区。这是因为对临界奇异性的严格处理需要采用非解析的(即非经典的)理论模型或方法, 如重整化群和跨接函数法, 但这些模型或方法极为复杂, 计算量也很大(郭涛等, 2016)。虽然一些半经验的方法可以近似地表达流体的临界奇异性, 但它们主要针对中低温压下的工业流体。由于上述原因, 加之N2的临界温度很低, 新方程的参数拟合与优化未考虑临界点附近的区域。

经过反复尝试, 新方程的参数拟合与优化能够确保方程精确地适用于273.15～4273.15 K。考虑到一些高温、高压研究的需要, 参数的优化力求让新方程能在高达5000 K左右仍然具有较高的精度。为了兼顾流体包裹体研究的需要, 参数的优化力求使所得方程能较好地外延到200 K附近。这一温度比CO2的三相点(216.592 K)还要低十几度。

新方程的参数拟合采用非线性最小二乘法。经过反复优化, 所得结果见表1。

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

3.3 对参考模型及其高温外延结果的再现

表1中的参数可以精确地再现参考模型的计算结果及其等容线的高温外延结果。根据再现结果, 方程(17)最适宜的条件是273.15～4273.15 K、0～0.85 g·cm−3。方程(17)在273.15～2000 K区间内的再现结果见图1, 在2000～4273.15 K区间内的高温外延的再现结果见图2。

图1中的最大体积偏差(the maximum deviation of molar volume, MDV)为0.306%, 平均体积偏差(average deviation of molar volume, ADV)为0.089%。总体上看, 在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; 参考模型数值据Span et al., 2000。

图2 方程(17)与参考模型等容外延结果的比较(=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以上的外延结果也与方程(17)吻合得很好(图3a), 其中MDV=0.60%, ADV= 0.27%。方程(17)向低温方向适当的外延结果如图3(b)所示, 其中MDV<2%。

由于立方型方程能够应用的最高密度max比较有限, 当密度接近或超过max时, 方程的偏差会随着密度的增大而迅速增大。当密度范围一定时, 方程的最大偏差(或偏差的变化范围)也是一定的。根据方程(17)的最大偏差, 本研究划分出几个密度区间(表2), 并按区间统计ADV与MDV, 其中也包括273.15 K以下和4273.15 K以上的偏差统计结果。

在273.15～4273.15 K区间内, 当max为0.850、0.926、0.964、1.000和1.030 g·cm−3时, 方程(17)相对于参考模型的MDV分别为0.41%、0.95%、1.39%、2.00%和2.50%。为了便于应用, 本文用方程(17)计算了这5个密度对应的等容线(图4), 并用经验式(21)对它们进行了拟合(表3)。

对于每条等容线, 压力随着温度的增大而增大。因此, 在方程能够应用的max以内, 温度越低, 方程能够应用的最大压力(max)也就越低。尽管如此, 方程(17)的适用范围基本上涵盖了一些典型的地温、地压梯度线(Tarney and Windley, 1977; Lambert, 1983; Chapman, 1986)及典型的变质岩相的稳定区(Dipietro, 2013), 包括高温、超高温变质作用的温压条件(刘守偈等, 2008; 魏春景, 2016; 李旭平等, 2019)。因此, 对于不同的地壳温压梯度, 方程(17)可适用于下地壳乃至上地幔的温压条件。

4 新方程的检验

4.1 与高温、高压实验数据的比较

现有的高温、高压实验数据并不多。本文选取了5组具有代表性的高温、高压实验数据(Tsiklis and Polyakov, 1968; Robertson and Babb, 1969; Malbrunot and Vodar, 1973; Mills et al., 1975; Antanovich and Plotnikov, 1976)与新方程进行了比较, 结果如图5和表4所示。其中Robertson and Babb (1969)的数据中约有1/3超过了方程(17)最适宜的密度范围(0～0.85 g·cm−3); Mills et al. (1975)的数据大部分都超过了0.85 g·cm−3, 所以在此只能选取密度在0.85 g·cm−3以下及附近(0.80～0.97 g·cm−3)的数据。这些数据相对于参考模型计算值的偏差见表4。其中Malbrunot and Vodar (1973)的体积数据在700 K以上有明显的负偏差, 并且随着温度的升高而增大, 最大值达−5%, 平均约为−2%。

从总体上看, 新方程在其适用范围内与实验数据均吻合得很好。其中Malbrunot and Vodar (1973)在800 K以上的体积数据由于有负偏差而均明显地低于方程(17)的计算值。与Mills et al. (1975)的数据相比, 方程(17)的偏差总体上较大。这是因为其多数数据的压力(或密度)均显著地超过了方程(17)最适宜的范围。不过, 随着温度的升高或压力的降低, 方程(17)的偏差下降得很明显。这与方程(17)对拟合数据的再现偏差的变化趋势是一致的。

cal代表方程(17); ref代表参考模型; 参考模型数据据Span et al., 2000。

表2 方程(17)的体积偏差与温度、密度范围的关系

注: ADV为平均体积偏差, MDV为最大体积偏差。

4.2 与高温、高压MD、MC模拟数据的比较

目前, 虽然原子–分子尺度上的计算机模拟提供了很多高温、高压下N2的数据, 但其中只有少数几组数据可与新方程做比较。例如, Etters et al. (1986)将的数值计算结果拟合为函数, 同时加上多极相互作用和分子内相互作用的贡献, 得到了N2的分子间势能函数, 并根据这一函数用恒压MC方法模拟了N2的性质(300～1500 K, 0.2～15 GPa)。Strąk and Krukowski (2007)采用高精度的量子力学从头算方法得到了N2的分子间势能, 并将这些势能用于N2的性质的MD模拟(300～2000 K, 0～40 GPa)。这两组数据都有很多高密度数据, 但本文只需要其中≤1.2 g·cm−3的数据, 其范围见表5。

Etters et al. (1986)数据中0.85 g·cm−3的数据点约占2/3, 其中最大密度达1.17 g·cm−3; Strąk and Krukowski (2007)的数据中0.85 g·cm−3的数据点也不少, 其中最大密度达到了1.20 g·cm−3。

上述两组模拟数据的最高压力为4.43 GPa, 而最高温度只有2000 K。为了检验方程(17)在高温、高压下的精度, 本研究应用Materials Studio软件的Forcite模块在更宽的-范围内(1473.15～5073.15 K, 0.3～6.0 GPa)对中高密度(0.4～1.196 g·cm−3)的N2进行了MD模拟。

表3 等容线方程(21)的参数

注: E+表示×10+n。

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

模拟采用系综(固定分子数的等温等压系综), 其对象为3000个N2分子, 其中N2的分子间势能取自COMPASS 库(Yang et al., 2000)。每个宏观状态点的模拟分成几个阶段, 微观状态之间的时间步长均为1 fs。第一阶段是初步模拟, 时长为50 ps。实际上系统趋于热力学平衡只需要20～40 ps。在后面的每个阶段, 时长都是1000 ps, 共计106个微观状态, 其中初态为前一阶段模拟的终态。每个阶段的模拟结束后,和的系综平均值及不确定度都会被记录。如果最近两个阶段的平均值之差在前一阶段的不确定度范围内, 则停止模拟。最后结果取最后一个阶段的平均值。在中高密度(>0.3 g∙cm−3)条件下, 前三个阶段的计算即足以达到模拟要求。模拟的范围见表5。比较结果表明, 这些数据与参考模型或其高温外延结果均吻合得很好。

从图6可以看出, 表5中的3组模拟结果与新方程均吻合得较好。其中较大的偏差都出现在高密度区, 相应的ADV和MDV也列于表5。

表4 一些高温、高压实验数据及其相对于参考模型的体积偏差(参考模型据Span et al., 2000)

注: ADV为平均体积偏差; MDV为最大体积偏差;为数据点的数量。

表5 方程(17)与一些高温、高压PVT或PρT模拟数据的比较

注: ADV为平均体积偏差; MDV为最大体积偏差;为数据点的数量。

图6 方程(17)与分子动力学、Monte Carlo模拟数据的比较(图a数据据Etters et al., 1986; 图b数据据Strąk and Krukowski, 2007)

5 与其他地质流体状态方程的比较

5.1 与参考模型及其高温外延结果的比较

HP91方程并没有说明其适用范围。此方程与参考模型及其高温外延结果的比较见图7。结果表明, 影响HP91方程精度的主要因素是温度, 这是因为参数与的关系过于简单。在400 K以上, 该方程的体积偏差随着温度的上升以近似线性的方式快速增大, 到1300 K附近达4%左右(图7a)。这些结果表明, HP91方程的应用不宜超过1300 K。在1300 K以下, 该方程与参考模型的比较见图7b。其中, 体积偏差随着压力的增加有明显的下降趋势。

据报道, DMW92方程可用到2000 K、2.0 GPa, 平均偏差约为1.5%。为了厘清DMW92方程实际的适用范围, 本研究将该方程与参考模型及其高温外延结果进行了比较(图8)。结果表明: ①压力是影响该方程精确度的主要因素, 其中压力与体积偏差之间具有显著的线性相关性; ②当压力增至0.7 GPa时, MDV约为3%左右, 而相应的温度可达4273.15 K左右。当压力增至1.0 GPa时, MDV增大到5%左右。

BS92方程据称可用于400～4000 K、0.5～100 GPa。该方程与参考模型及其高温外延结果的比较如图9所示。很明显, 在0.8 GPa以下, 体积偏差随着压力的下降迅速增大, 其中MDV达9.54%; 在1.2 GPa以上, 体积偏差随着压力的增大也上升较快, 在2.9 GPa附近达到6.40%。

据原文报道, B99方程可应用到1300 K和1 GPa。本文用参考模型及其高温外延结果对它进行了验证。图10可见, 体积偏差随着压力的上升以近似于直线的方式快速增大, 在0.75 GPa和1 GPa时达到5.01%和7.34%。此方程可以应用的温度远高于1300 K,但必须限制压力的范围。

5.2 与高温、高压MD、MC模拟数据的比较

HP91、DMW92、BS92与B99方程与表5中的MD、MC模拟数据的比较结果见图11。比较时对各方程及模拟数据的限制条件如下, HP91:≤1300 K,≤2.5 GPa; BS92:≥0.5 GPa; B99:≤1.0 GPa; DMW92:≤2.0 GPa; 方程(17)没有限制。这些条件要么在各方程原文声明的范围内, 要么在本研究所确认的适用范围内。

图7 HP91方程与参考模型的对比(cal数据据Holland and Powell, 1991; ref数据据Span et al., 2000)

图8 DMW92方程与参考模型及其高温外延结果的比较(cal数据据Duan et al., 1992; ref数据据Span et al., 2000)

图9 BS92方程与参考模型及其高温外延结果的比较(cal数据据Belonoshko and Saxena, 1992; ref数据据Span et al., 2000)

图10 B99方程与参考模型及其高温外延结果的比较(cal数据据Bakker, 1999; ref数据据Span et al., 2000)

从总体上看, HP91、DMW92与B99方程的偏差均较大或很大, 方程(17)的偏差最小, 其次是BS92方程。与Etters et al. (1986)、Strąk and Krukowski (2007)和本研究的模拟结果相比, 方程(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)的相应偏差, 其主要原因是BS92方程在1 GPa以下的偏差较大。

图11进一步印证了方程(17)及上述4个方程与参考模型及其高温外延结果的比较结果。这也表明, 方程(17)可以扩展到比图4中的等容线更高的条件, 但必然伴随着偏差的上升。

根据方程(17)计算的等容线, 当≤1.2 g·cm–3时, 方程(17)在200～5573.15 K范围内的结果都是合理的。结合本研究的结果可以推测, 当≤1.2 g·cm–3时, 如果≤5073.15 K, 方程(17)的MDV约为4%, 相应的压力不超过6.8 GPa; 如果≤5573.15 K, 方程(17)的MDV约为4.5%, 相应的压力不超过7.5 GPa。

以上结果表明, 方程(17)在其适用范围内均明显优于HP91、DMW92、BS92与B99方程。其中DMW92与B99方程的适用范围与原文报道的结果相差很大, 而BS92方程在0.8 GPa以下偏差较大。

6 其他热力学性质的预测

含氮物质的迁移、交换、转化过程(如扩散、不同相间的物质交换、相变、化学反应)取决于含氮物质的各种热力学与动力学性质等。这些过程也伴随着能量的交换、转化和力的传递。要想正确地理解有流体参与的地质过程, 不但需要计算流体的性质, 还需要计算流体的膨胀系数()、压缩系数()、声速()、等压或等容摩尔热容(C,C)、焓()、熵()、Gibbs自由能()、Helmholtz自由能()、内能()、化学势()、逸度()、活度()等许多热力学性质。例如, 相平衡、化学平衡、热效应的计算则需要用到、、、、、等。另外, 一些流体动力学性质的计算也涉及热力学性质。例如, 流体的声速和粘度()的计算需要用到流体的等压及等容热容、密度、膨胀系数和压缩系数等。

上述的各种热力学性质都可以根据状态方程、热容、标准焓、标准熵以及相关的热力学关系式推算出来。根据状态方程, 很容易计算出膨胀系数、压缩系数、逸度系数()、摩尔剩余焓(R)、摩尔剩余熵(R)等性质。在此基础上可以进一步推算化学势、逸度、活度等许多其他的热力学性质。实际上,、、R、R、或等性质很少有实验或模拟数据的报道。不过, 通过参考模型可以计算这些性质。在此基础上, 根据剩余性质的定义及理想气体的焓、熵与压力的关系, 用参考模型算出真实流体的R和R, 然后很容易算出R和ln。

根据附录A中的表达式和方程(17)预测ln、R、R和C–C。其中应该采用方程(17)在给定温压下的体积根(附录B)。这些性质的预测结果与参考模型计算值的比较见图12。可以看出, 方程(17)与参考模型对ln和R的计算结果吻合得非常好。即使在低温区(200～273.15 K), ln的计算结果也吻合得很好。R和C–C的计算结果大多数吻合得很好。即使C–C在275 K附近的几个峰也被精确地重现。实际上, 由于状态方程推导C–C的过程涉及与对温度的二阶导数, 所以C–C对状态方程的偏差很敏感, 其计算精度可以作为检验状态方程精度的一个很好的指标。

(a)sim为MC模拟数据, 据Etters et al. 1986; (b)sim为MD模拟数据, 据Strąk and Krukowski, 2007; (c)sim为本研究的MD模拟数据。

图12 方程(17)的预测结果与的参考模型计算值的比较(200～2000 K)(参考模型据Span et al., 2000)

不过, 在低温区, 随着压力的增大,R和C–C的偏差逐渐增大。这主要是因为方程(17)的偏差随着温度的降低而逐渐增大。

此外,R在低温、高压区的偏差较大还有另外的因素。设Δln、ΔR与ΔR分别为ln、R与R计算结果的偏差。由于R=R–R,R=ln, 从而可得ΔR–ΔR=Δln。另一方面, 由于ln的计算结果很精确, 可以假设Δln≈0。由此可得ΔR≈ΔR/。因此, 当ΔR相同时, 温度越低, ΔR越大。也就是说, 低温放大了R的偏差。不过, 这几乎不影响Gibbs自由能计算结果的精确性。因为ln的计算结果很精确, 所以R与R的偏差几乎总能相互抵消。

7 结 论

为了精确地计算高温、高压下含N2地质流体混合物的热力学性质, 本研究根据N2的高精度数据建立了超临界氮的一种高精度四参数立方型状态方程。在273.15～4273.15 K、0～0.85 g·cm−3范围内, 该方程的体积偏差在0.41%以内, 平均体积偏差不到0.1%。新方程可以比较精确地应用到5073.15 K或0.926 g·cm−3, 相应的体积偏差在0.6%或1.0%以内。新方程可以较好地外延至200 K、1.03 g·cm−3, 相应的体积偏差不超过2.5%。新方程也可近似地应用到1.2 g·cm−3, 相应的最高温压可达5573.15 K、7.5 GPa。

为了检验新方程在高温、高压下的精度, 本文利用分子动力学方法模拟了N2在1473.15～5073.15 K、0.3～6.0 GPa范围内的性质。这些结果及现有的高温或高压实验数据、分子动力学及Monte Carlo模拟数据(0.1～1.2 g·cm−3)均与新方程吻合得很好或较好, 相应的偏差小于或接近于实验或模拟的不确定度。由新方程预测的膨胀系数、压缩系数、逸度系数、剩余焓、剩余熵等热力学性质均与参考模型的计算值吻合得很好。这些结果表明, 新方程显著地优于现有的其他地质流体状态方程, 从而为高温、高压下多元含氮混合物热力学性质的精确模拟或预测提供了一个较好的纯组分模型。

新方程简单、易用, 外延性也较好, 也易于向混合物推广, 但在精确度上不可能超过参考模型。另外, 参考模型参数拟合所用的部分高温或高压实验数据存在一定的系统偏差或较大的不确定度, 其可靠性明显不如中低温压下的实验数据。未来的高温、高压实验若能提供新的高精度数据, 则新方程的参数还可以进一步优化。

附录A 一些基本热力学性质的表达式

A.1 膨胀系数与压缩系数

由于只是的隐函数, (A1)式中第二个等式的变换采用了隐函数的求导法则。考虑到文献中一般都没有与数据, 本研究采用C–C来间接地检验与计算结果的准确性:

A.2 逸度系数

在计算流体组分的化学势(即摩尔Gibbs自由能)之前, 必须先计算组分的逸度系数:

在方程(17)中, 由于是和的显函数, 采用公式(A11)来推导ln更为方便:

式中=/(压缩因子)。逸度系数及剩余焓、剩余熵的表达式需要对状态方程进行积分才能得到。积分结果的形式取决于方程(17)中二次项分母(+)的根的情况。(+)的根的判别式Δ=2。由于Δ可能因温度的变化而改变符号, 所以Δ的值有两种情况: 当≠0时, Δ>0; 当=0时, Δ=0。当Δ>0时, 将方程(17)代入公式(A11)可得:

当Δ=0时,

A.3 摩尔剩余焓与摩尔剩余熵

对于显压型方程[=(,)]可以直接由状态方程推导摩尔剩余内能R(即摩尔剩余热力学能)。然后再由公式(A15)算出摩尔剩余焓R:

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

当Δ>0时, 将方程(17)代入公式(A16), 积分后可得

当Δ=0时, 将方程(17)代入公式(A16), 积分后即可得到R的表达式。当Δ=0时,=0。据此可将R的表达式简化为

当Δ=0时, 温度只是一个特定的值(温度轴上的一个点)。因此可以说, 出现Δ=0的概率无穷小。在具体的计算中实际上无需考虑这种情况。

附录B 方程(17)在给定-条件下的求解

首先考虑解析法。任一立方型状态方程都可统一转化成如下形式:

其中可以是压缩因子()、摩尔体积()或数密度(即的倒数)。然后, 将系数1～3代入求根公式即可得到体积根, 详见郭涛等(2016)的附录A。

其次, 也可以采用二分法或其他方法。其中二分法是求解各种非线性方程最常用的方法之一。对于一定温度下的超临界流体, 由状态方程(17)算出的关系是一条单调的曲线。因此, 方程(17)在其适用范围内的任意-条件下只有唯一的体积根, 满足二分法的应用条件。

在给定-条件下, 应用二分法求体积根的思路大体如下:

①首先定义一个误差函数(,), 如(,)=cal(,)/1。其中cal是由方程(17)算出的压力,为给定的压力,是一个任意的摩尔体积。与此同时, 为迭代计算设定一个允许的最大误差, 如=1.0×10−8。

②为待求根设置两个合适的初值(即1和2), 以便粗略地圈定待求根所在的区间[1,2]。如果取初值时没有任何参考值可用, 可取1=,2=(1+)。其中为方程(17)的协体积,为一个极小的正值;可取0～1之间某个较小的正值, 如0.1或0.2。如果有实验值或其他值可作为参考(ref), 可令1=(1–)ref,2=(1–)ref, 其中ref为参考值。

③将温度与初值1和2代入方程(17), 得到1=cal(,1),2=cal(,2), 然后算出两个压力下的偏差, 即1=(,1)＝1/1,2=(,2)=2/−1。

④若1·2>0, 说明方程(17)在[1,2]内无根, 然后返回步骤②, 调整初值, 使其满足1·2<0; 若1·2<0, 说明方程(17)在[1,2]内有根, 然后进入下一步的迭代循环。

⑤令3=(1+2)/2, 算出相应的压力3=cal(,3); 同时算出相应的压力偏差3=(,3)=3/−1。如果|3|<, 则3就是允许误差范围内待求根的数值解。此时可直接跳到步骤⑧; 否则, 执行下一步运算。

⑥如果|3|>, 则更新[1,2]的边界。具体方法: 若1·3>0, 即=1和=3时的压力偏差同号。这说明3比1更接近于待求根。此时应将边界值1更新, 即令1=3,1=3,1=3。若1·3<0, 则进行相反的运算, 即令1=2,1=2,1=2。

⑦重复步骤⑤和⑥。

⑧记录迭代结果, 结束求根。

致谢:本文作者对清华大学化学工程系于养信教授在分子动力学模拟方面的指导和大力支持深表感谢, 同时感谢两位匿名审稿专家提出的建设性修改意见。

段振豪. 2010. 地质流体状态方程. 中国科学: 地球科学, 40(4): 393–413.

郭磊. 2017. 高温高压下二氧化碳状态方程的评价. 石家庄: 河北地质大学硕士学位论文: 1–50.

郭涛, 胡家文, 王向辉, 靳丽花, 王艳哲. 2016. 超临界水的热力学模拟: 一种可用到4273 K、2 GPa的立方型状态方程. 岩石学报, 32(7): 2196–2208.

胡家文. 2002. 维里型和立方型地质流体状态方程的理论研究. 成都: 成都理工大学博士学位论文: 1–101.

黄荣荣. 1995. Burnett法--测定装置. 江苏石油化工学院学报, 7(3): 8–13.

李旭平, 王晗, 孔凡梅. 2019. 高温–超高温变质作用成因研究——来自华北克拉通西部孔兹岩带和南非Kaapvaal克拉通西南部Namaqua活动带与Bushveld变质杂岩体的启示. 岩石学报, 35(2): 295–311.

刘斌. 2011. 简单体系水溶液包裹体pH和Eh的计算. 岩石学报, 27(5): 1533–1542.

刘守偈, 李江海, Santosh M. 2008. 内蒙古土贵乌拉孔兹岩带超高温变质作用: 变质反应结构及指示. 岩石学报, 24(6): 1185–1192.

卢焕章, 范宏瑞, 倪培, 欧光习, 沈昆, 张文淮. 2004. 流体包裹体. 北京: 科学出版社: 1–485.

倪培, 迟哲, 潘君屹, 王国光, 陈辉, 丁俊英. 2018. 热液矿床的成矿流体与成矿机制——以中国若干典型矿床为例. 矿物岩石地球化学通报, 37(3): 370–394.

孙贺, 肖益林. 2009. 流体包裹体研究: 进展、地质应用及展望. 地球科学进展, 24(10): 1105–1120.

王艳哲, 2015. 高温高压下水的状态方程的检验与评价. 石家庄: 石家庄经济学院硕士学位论文: 1–71.

魏春景. 2016. 麻粒岩相变质作用与花岗岩成因–Ⅱ: 变质泥质岩高温–超高温变质相平衡与S型花岗岩成因的定量模拟. 岩石学报, 32(6): 1625–1643.

肖益林, 余成龙, 王洋洋, 陆一敢. 2018. 变质作用与流体包裹体: 进展与展望. 矿物岩石地球化学通报, 37(3): 424–440.

徐航, 蔡旭东, 胡芃, 金熠. 2015. 基于磁悬浮密度计的高精度流体测量系统. 中国科学技术大学学报, 45(5): 404–408.

杨光璐, 蔺玉秋, 刘加林, 鲍君刚. 2004. 辽河油田稠油油藏氮气泡沫驱适应性研究. 新疆石油地质, 22(2): 188–190.

翟伟, 孙晓明, 徐莉, 张泽明, 梁金龙, 梁业恒, 沈昆. 2005. 苏北青龙山超高压变质榴辉岩流体包裹体特征与流体演化. 岩石学报, 21(2): 482–488.

张泽明, 沈昆, 赵旭东, 石超. 2006. 超高压变质作用过程中的流体——来自苏鲁超高压变质岩岩石学、氧同位素和流体包裹体研究的限定. 岩石学报, 22(7): 1985–1998.

周健, 陆小华, 王延儒, 时钧. 1999. 超临界水的分子动力学模拟. 物理化学学报, 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 (Cristalândia, 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

河北省自然科学基金(D2018403089)、河北省高等学校科学技术研究重点项目(ZD2018026)和河北省研究生创新资助项目(CXZZSS2020112)联合资助。

姜思雨(1996–), 男, 硕士研究生, 地球化学专业。E-mail: jiang_geol@sina.com

胡家文(1967–), 男, 教授, 主要从事计算地球化学研究。E-mail: hu_jiawen@sina.com