超临界态二氧化碳流体热力学性质的分子动力学模拟

刘丽丽 代威 杨智 廖传军

(1 西安航天神舟建筑设计院有限公司北京分公司 北京 100040)

(2 北京宇航系统工程研究所 北京 100076)

(3 广东工业大学材料与能源学院 广州 510006)

1 引言

二氧化碳化学性质稳定,无毒,不易燃易爆,臭氧消耗和气候变暖潜值较低,天然绿色环保,临界状态易于实现。在临界点以上,二氧化碳具有类似于液体的密度,气体的粘度和扩散系数,溶解能力强,表面张力为零,具有良好的传热传质特性[1]。目前,CO2在超临界萃取、制冷热泵及热力发电系统、页岩压裂及强化石油开采和环境保护(如CO2捕集和封存)等诸多领域已有广泛应用。然而,超临界态流体热物性的准确获取是其得以广泛应用的重要前提。

传统上,流体热物性可通过实验测量以及状态方程的关联计算来获取。然而,为获取超临界态热物性数据,实验测量方法则需针对测量装置的耐高温高压性能进行特殊设计,同时还需保证较高的测量精度并测量足够多的数据点,其代价通常极为昂贵且技术上难以实现。对状态方程而言,其通常对实验数据依赖性较强,如专用型状态方程[2]的确定则需较多的参数自由度及大量的实验数据样本,其适用范围也严格受限于实验数据样本所处的热力状态区间,外推预测能力较弱,而简单的立方型状态方程通常外推预测精度相对不足。

本质上,流体的宏观热物性由流体分子的微观相互作用(分子力场)所决定。作为热物性获取的一种有效方法,基于力场的分子模拟具有坚实的统计热力学理论基础。近年来,随着计算机技术的快速发展,分子模拟在流体热物性研究领域扮演越来越重要的角色,其不仅可提供分子体系微观结构的可视化研究,实现流体宏观热物性的微观机理解释,也可作为实验测量方法的一种有效补充,获得具有定量精度的内插及外推物性数据、获得实验无法测量或难以发现的结果(如极端的高温高压条件)、预测现有或新材料的性质及发展新理论。分子模拟在流体热物性研究领域已被公认为联系实验数据和理论模型的重要桥梁。

分子力场是分子模拟的关键,目前国内外学者针对CO2流体已开发出多种不同的分子力场,如简单的粗粒化单粒子模型[3]以及高分辨率的全原子模型[4-6]均得到了广泛应用。尽管CO2分子力场存在诸多不同的描述形式,但大多数力场中描述分子间范德华相互作用的参数均通过拟合相平衡实验数据来获取,以实现其对亚临界态热物性较好的预测精度。然而据文献调研发现,基于气液相平衡数据拟合的半经验力场模型对超临界态热物性的外延预测性能仍少有研究。因此,采用分子动力学模拟方法,针对几种典型的粗粒化和全原子力场模型,如SAFT-γ[3]、Zhang[4]、EMP2[5]、TraPPE[6]等,开展超临界态CO2流体热力学特性的预测研究,以充分评估不同力场对CO2流体超临界态热力学性质的外延预测性能,从而为极端温度和压力工况下CO2流体热物性的获取提供一种有效的替代方案。

2 模拟方法

2.1 力场模型

本研究中,不同力场的范德华相互作用均可采用一般化Mie 势能函数[7]进行描述,其函数形式为:

对于带电粒子模型,分子间还需考虑长程静电相互作用,如式(3),其中qi和qj为位于原子中心的点电荷,ε0为真空介电常数。对于全原子力场模型,体系势能还包含分子内原子间相互作用,也即键伸缩项和键角弯曲项,如式(4)。其中,kr、kθ分别为键伸缩和键角弯曲力常数,r0、θ0分别为平衡键长和键角。主要考虑了4 种CO2分子力场,其中TraPPE、EPM2和Zhang 模型为三点硬球模型,SAFT-γ 模型为单点粗粒化模型。相比于三点模型,单点模型通常具有更简单的结构和更高的计算效率,在相同的计算条件下,可用于更大尺度流体行为和现象的模拟。4 种CO2力场参数如表1 所示;TraPPE 和Zhang 模型采用Lorentz-Berthelot 混合规则描述不同类原子间相互作用,而EPM2 力场采用几何平均数规则。

表1 不同CO2 力场模型对应的相互作用参数Table 1 Interaction parameters of different CO2 field models

2.2 模拟细节

所有模拟均在开源的LAMMPS 平台下进行,对于单点粗粒化力场模型,模拟体系包含1 000 个粒子,对于全原子力场模型,模拟体系由512 分子填充。模拟均在恒温恒压(NPT)系综下进行,且采用周期性边界条件。模拟的时间步长设置为1 飞秒,同时采用Nosé-Hoover 方法对模拟体系进行温度和压力控制。在数据采集阶段,采用PPPM 方法计算长程库伦相互作用,且将静电力计算的相对偏差控制在1 ×10-4以内。每个工况点下,模拟的总步数设置为750 万步,其中前500 万步用于系统的弛豫平衡,后250 万步用于数据采用及热力学性质的统计平均。

2.3 热物性参数计算

在二氧化碳的热物性参数计算中,流体的密度ρ可直接通过系统体积的统计平均值来计算,也即:

式中:N为模拟体系的分子数,M为分子摩尔质量,V为模拟系统体积,Na为阿伏伽德罗常数。

通常,流体的比定压热容Cp由两部分组成,也即理想项和剩余项:

理想项贡献Cid p取自于实验关联式,剩余项贡献则是通过系综平均获得,其表达式如式(7)所示:

式中:KB为玻尔兹曼常数,T和P分别为体系所设定的温度和压力,Hconf、Uconf分别为体系的构型焓和构型内能(也即排除分子动能和分子内原子间的相互作用)。

3 结果与讨论

3.1 密度(ρ)预测

采用4 种典型的分子力场模型(TraPPE、EPM2、Zhang、SAFT-γ Mie),采用分子动力学模拟方法,在NPT 系综下,对CO2流体在超临界态工况条件下(T=600—900 K,p=30—100 MPa)的密度特性进行了模拟预测,且将模拟预测结果与NIST 数据库中的参考数据进行了对比,如图1 所示。同时,不同温度下各分子力场密度预测的平均绝对相对误差如表2所示。由图1 和表2 可知,4 种力场模型对高温高压条件下CO2流体的密度预测值与NIST 参考值的平均绝对相对偏差均在3%以内,因此均具有较好预测精度。值得注意的是,3 种全原子力场(TraPPE、EPM2、Zhang)对密度的预测精度均随着温度的升高呈现逐渐减小的趋势,也即在高温条件下具有更好的预测精度,而单点粗粒化SAFT-γ Mie 力场密度预测精度随温度的变化趋势却相反。尽管如此,简化的SAFT-γ Mie 力场总体上在超临界态工况条件下(T=600—900 K,p=30—100 MPa)对密度具有相对最好的预测精度,这主要由于SAFT-γ Mie 力场模型参数的确定采用了一种全局化的参数优化策略。

图1 不同力场所预测的CO2 流体密度值与NIST 参数数据的对比Fig.1 Comparison of predicted CO2 fluid density values using different force fields and NIST parameter data

表2 不同力场在不同温度下密度预测的平均绝对相对偏差Table 2 Average absolute relative deviation of density prediction using different force fields at different temperatures

3.2 比定压热容(Cp)预测

同样在NPT 系综下,对超临界态工况下(T=600—900 K,p=30—100 MPa)CO2流体的比定压热容特性进行了模拟预测,模拟值与NIST 参考数据对比结果如图2 所示。其中,表3 为不同温度下各分子力场所预测的平均绝对相对误差。由图2 和表3 可知,针对超临界态比定压热容,4 种力场模型所预测的平均绝对相对偏差均在1%以内,因此比密度的预测精度更高。由表3 可知,EPM2 和Zhang 力场模型对比定压热容的预测偏差均随温度的升高而逐渐减小,而TraPPE 和SAFT-γ Mie 力场的预测精度随温度的变化不显著。总体上,在温度T=600—900 K 和压力p=30—100 MPa 的超临界态范围内,EPM2 力场对CO2流体的比定压热容特性具有最好的预测精度。

表3 超临界工况下,不同力场在不同温度下的比定压热容预测偏差Table 3 Average absolute relative deviation of specific heat capacity at constant pressure using different force fields at different temperatures

图2 不同力场所预测的CO2 流体比定压热容与NIST 参数数据的对比Fig.2 Comparison of predicted CO2 specific heat capacity at constant pressure using different force fields and NIST parameter data

同时,上述模拟结果表明,基于气液相平衡特性拟合关联的CO2力场模型对超临界态热力学性质(如密度和比定压热容)同样具有较好的预测精度。此外,尽管全原子力场模型具有更多的参数自由度,但在超临界态热力学性质预测中,其预测性能相比粗粒化力场模型并非具有显著优势。由于在温度T=600—900 K 和压力p=30—100 MPa 的超临界态区间内,CO2流体不具有显著的变物性行为,因此4 种力场模型均体现出较好的预测精度。同时,随着温度的升高,EPM2 和Zhang 两种力场对密度和定压比热容的预测偏差均呈现出逐渐减小的趋势,这表明该力场的预测能力可推广应用于更高温度的热力工况。

4 结论

采用分子动力学模拟方法结合4 种典型的半经验型分子力场,对超临界态二氧化碳流体热力学性质进行了模拟预测研究,并得出以下结论:

(1)尽管4 种力场的范德华作用参数通过拟合气液相平衡数据而获取,但在温度T=600—900 K 和压力p=30—100 MPa 的超临界态区间内,4 种力场对流体密度和比定压热容特性均具有较好的预测精度,这表明4 种力场模型均具有较好的外延预测性能。

(2)EPM2 和Zhang 力场在对密度和比定压热容的模拟预测中,其预测偏差随温度升高均呈现出逐渐减下的趋势,这表明该力场可用于更高温度工况条件下CO2流体热力学性质的预测。

(3)相比于简单且计算高效的单粒子SAFT-γ Mie 力场,全原子力场模型尽管具有更多的参数自由度,但其在超临界态区间内对CO2流体热力学特性的预测并不具有显著优势。这表明全原子力场模型的力场参数对于超临界态热物性预测而言并非最优解。