超临界CO2体系黏度分子动力学研究

靳遵龙， 李德雨， 刘东来, 王永庆， 王定标

(1.郑州大学 化工与能源学院 河南 郑州 450001; 2.洛阳双瑞特种装备有限公司 河南 洛阳 471003)



超临界CO2体系黏度分子动力学研究

靳遵龙1,2， 李德雨2， 刘东来1, 王永庆1， 王定标1

(1.郑州大学 化工与能源学院 河南 郑州 450001; 2.洛阳双瑞特种装备有限公司 河南 洛阳 471003)

由于具有较低的临界温度以及适中的临界压力，超临界CO2成为最为常见的超临界流体，被广泛应用到超临界萃取和制冷等领域.充实超临界CO2体系的基础物性数据，对于进一步利用超临界CO2具有积极意义.运用平衡分子动力学方法，以Lammps免费开放源代码模拟平台为基本框架，自行编译针对超临界CO2特性的脚本文件，对超临界CO2体系的黏度物性进行计算，并将计算结果与成熟的物性数据库进行对比.结果表明，所编译脚本文件计算值与相关数据误差在0.5%以内，说明使用该方法用来获得超临界CO2黏度是可行的.

超临界CO2； 分子动力学； 黏度

0 引言

超临界流体技术与应用受到了国内外众多学者和工业界的高度关注.超临界流体既不是气体又异于液体，具有许多独特的物理化学性质，其工程背景涵盖多个工业领域[1].

超临界CO2流体具有绿色环保等突出优点，被广泛应用到超临界萃取和制冷等领域[2].传统的制冷剂对臭氧层具有强烈的破坏性，且会引发温室效应.新型制冷剂具有不破坏臭氧层、几乎不引发温室效应、无毒无害、不易燃易爆、经济安全、资源丰富等特点.CO2作为天然制冷工质，具备上述特点的同时，兼备其他优点：优良的热物性及流动性、较高的单位容积制冷量、蒸发潜力大等，被广泛应用到众多工业领域.由于超临界CO2分子间相互作用的复杂性,现有的应用于气、液相的理论模型,在预测超临界CO2的性质时存在较大的难度.高温高压的实验条件也限制了超临界流体的实验研究.因此，超临界CO2体系的基础数据相对缺乏，这阻碍了超临界CO2的广泛应用.分子动力学模拟不受实验条件的限制,并随着计算机软、硬件的快速发展，对超临界CO2的物理化学性质进行分子动力学模拟研究,近年来得到了越来越广泛的应用[3].

文[4]应用分子动力学手段研究了温度从313～333 K，压力15 MPa下超临界CO2的无限稀释扩散系数，所得结果与实验结果吻合较好，表明分子动力学手段是一种用来获取超临界流体的基础物性数据的强大方法.文[5]采用从头算分子动力学模拟手段研究了丙酮-CO2混合体系，考察了在有限温度下弱电子给体-受体分子间相互作用的动力学影响.计算结果验证了静态电子计算方法所做的某些假设.文[6]结合Lennard-Jones流体对比状态方程, 按照对应态原理得到38种有机物分子的新的Lennard-Jones位能参数，并将这些参数应用于若干种有机物在超临界CO2中的无限稀释系数的分子动力学模拟当中，结果表明，利用这种思路可以有效地预测超临界体系的扩散性质.文[7]运用Car-Parrinello分子动力学方法研究了在超临界CO2中形成的邻羟基苯甲酸-水分子的结构.研究了原子间的径向分布函数、配位数以及平均氢键数等.文[8]运用分子动力学模拟方法研究了准临界、超临界CO2中正构烷烃的无限稀释系数, 模拟计算结果很好地与实验值吻合.文[9]运用分子动力学方法，研究了黏土蒙脱石在CO2插层存在情况下的聚合物改性.研究结果表明，散装聚合物能够混合吸收大量的超临界CO2.文献[10]运用经典的分子动力学理论模拟计算了超临界CO2以及其与极性溶剂：水、甲醇、乙醇(摩尔分数为0.125，T=318 K，ρ=0.7 g/cm3)的混合体系的特性，研究揭示了超临界CO2的分子结构特性以及CO2与极性溶剂相互作用机理.

作者应用分子动力学模拟软件Lammps自行编译程序，计算超临界CO2的黏度并与相关数据比较，一方面验证程序可靠性与准确性，另一方面为后期研究超临界CO2体系输运特性提供基础数据，为超临界CO2的研究提供分子水平的依据和参考.

1 模拟方法及模拟细节

通常情况下，分子模拟有两种方法来获得输运系数数据，一种是平衡分子动力学方法(EMD)，它通过平衡时间相关函数的Einstein 关系式或Green-Kubo 关系式来计算输运系数[11-12]；另一种方法是非平衡动力学方法(NEMD)，这种方法是通过对所研究系统施加外扰动，依靠扰动所产生的响应来求解相关特性[13].本文应用平衡分子动力学方法，应用Lammps 模拟软件自行编译程序，计算超临界CO2的黏度.

Lammps 软件是由Sandia国家实验室开发的开源代码，可以对气体、液体和固体以及它们混合体系的原子、聚合物等，从几个粒子乃至上亿个粒子同时进行建模计算，它既可以单机运行，又可以并行机运算，后者运算特点更突出.Lammps软件是免费源代码，因此可以根据需要对其进行修改和扩展.

图1 CO2元胞结构图

利用Lammps进行计算需要一个脚本文件(input_script)，脚本文件的命名格式为*.in.运行时需要将这个文件与Lammps.exe 放到同一个文件夹下.

一个脚本文件需要4部分组成：

1) 设定初始值和初始状态(单位、维数、边界条件、粒子类型等)；

2) 粒子的定义(粒子空间坐标、盒子大小等)；

3) 设置(仿真参数、电场系数、输出选项等)；

4) 运行脚本(通过Run 命令实现).

具体计算时，应用Material Studio建立模型如图1所示，获得拓扑信息后生成*.data文件格式，导入Lammps，编译程序进行计算.计算超临界CO2输运性质采用Green-Kubo[11]公式，最后结构提取采用print格式输出.

应用Green-Kubo[11]公式来计算剪切黏度系数，剪切黏度系数与张量非对角元素的时间自相关函数之间的关系式为

(1)

应力张量计算式：

(2)

2 结果与讨论

计算512个超临界CO2分子不同温度与不同压力下的黏度，并将计算结果与REFPROP 软件计算结果进行对比，REFPROP 物性软件计算黏度值和真实值误差在0.5%范围内，因此对比结果具有一定的参考价值.

图2为9 MPa不同温度下CO2黏度变化曲线，从图中可以清楚看到，程序计算出来的黏度值比REFPROP物性软件计算值大，但整体变化趋势一致，即在准临界温度附近，黏度曲线斜率陡，经过该区域后，黏度变化不剧烈.图3 为320 K不同压力下CO2黏度变化曲线，从图中可以看出，压力在15～55 MPa范围内，程序计算的结果与物性软件计算的结果误差较小，而准临界压力附近以及高压力条件下，计算的误差相对较大.这是由于准临界温度附近，物性变化剧烈，程序所使用的势能函数以及截断半径所带来的误差所致，而高压部分，猜测可能是CO2体系分子键以及分子物性改变所致，但整体变化趋势一致.因此，应用自行编译程序的方法可以对超临界CO2的输运特性进行计算.

图2 9 MPa下CO2黏度随温度变化曲线

图3 320 K下CO2黏度随压力变化曲线

3 结论

本文应用平衡分子动力学方法，应用Lammps开源模拟器平台，自行编译针对超临界CO2特性的脚本文件，计算超临界CO2体系不同温度与不同压力下的黏度.所编译程序的计算结果与成熟的物性数据库REFPROP进行对比，误差在0.5%范围内.结果表明，应用自行编译程序的方法对超临界CO2的输运特性进行计算是可行的.这些“机器实验数据”对超临界CO2的广泛应用可以提供必要的基础数据.

[1] 王飞，张光华，王三军，等．Fe团簇在Fe(110)表面上扩散和结构稳定性的分子动力学研究[J]．郑州大学学报：理学版，2012，44(3)：74-80．

[2] 杨凤叶，刘敏珊，刘彤，等．微细管内超临界CO2对流换热研究[J]．郑州大学学报：工学版，2014，35(2)：96-99．

[3] 何岩. 超临界CO2体系的扩散性质和微观结构的分子动力学研究[D]. 天津：天津大学，2007.

[4] Wang Jinyang, Zhong Haimin, Feng Huajie. Molecular dynamics simulation of diffusion coefficients and structural properties of some alkylbenzenes in supercritical carbon dioxide at infinite dilution[J]. Journal of Chemical Physics, 2014, 140(10): 104501.

[5] Munoz-Losa A, Martins-Costa M T C, Ingrosso F. Correlated ab initio molecular dynamics simulations of the acetone-carbon dioxide complex: implications for solubility in supercritical CO2[J]. Molecular Simulation, 2014, 40(1): 154-159.

[6] 周健，陆小华，王延儒，等. 有机物在超临界CO2中无限稀释扩散系数的分子动力学模拟[J]. 化工学报，1999，50(4)：491-499.

[7] Petrenko V E, Antipova M L, Gurina D L. The structure of hydrated complexes of o-hydroxybenzoic acid in water-modified supercritical carbon dioxide: the Car-Parrinello molecular dynamics simulation[J]. Journal of Supercritical Fluids, 2014, 85: 1-5.

[8] Feng Huajie, Gao Wei, Sun Zhenfan. Molecular dynamics simulation of diffusion and structure of some n-Alkanes in near critical and supercritical carbon dioxide at infinite dilution[J]. Journal of Physical Chemistry，2013，117(41)：12525-12534.

[9] Krishnan M, Saharay M, Kirkpatrick R J. Molecular dynamics modeling of CO2and Poly(ethylene glycol) in montmorillonite: the structure of clay-polymer composites and the incorporation of CO2[J]. Journal of Physical Chemistry, 2013, 117 (40): 20592-20609.

[10]Gurina D L, Antipova M L, Petrenko V E. Structural features of binary mixtures of supercritical CO2with polar entrainers by molecular dynamics simulation[J]. Russian Journal of Physical Chemistry, 2013, 87 (10): 1662-1667.

[11]Castai G，Ford J，Vivaldi F,et al. One dimensional classical many-body system having a normal thermal conductivity[J]. Physical Review Letters，1984, 52(21): 1861-1864.

[12]Poetzsh R H，Better H. Interplay of disorder and anharmonicity in heat conduction molecular-dynamics study[J]. Physical Review B，1994, 50(21):15757-15763.

[13]Allen M P，Tlldesley D J. Computer Simulation of Liquids[M]. Oxford：Clarendon Press，1987.

(责任编辑：王浩毅)

Molecular Dynamics Simulation of Viscosity in Supercritical CO2

JIN Zun-long1,2, LI De-yu2, LIU Dong-lai1, WANG Yong-qing1, WANG Ding-biao1

(1.SchoolofChemicalEngineeringandEnergy,ZhengzhouUniversity,Zhengzhou450001,China;2.LuoyangSunruiSpecialEquipmentCO.,LTD.,Luoyang471003,China)

Supercritical CO2was the most commonly used supercritical fluid due to its low critical emperature and moderate critical pressure. It has found wide applications in various fields, such as supercritical fluid extraction, cooling system, and rapid expansion of supercritical solution, etc. Enriching the basic physic data of supercritical CO2is of positive significance. Equilibrium molecular dynamics simulations were carried out in this study to investigate the viscosity of supercritical CO2. The numeration process was conducted by our compiled script, based on the Lammps simulator, which was a classical molecular dynamics code. The calculation outcome were compared with the mature property database. The results showed that the error was less than 0.5% between the calculation and that of the database. And it was feasible to use the presented method to obtain the supercritical CO2viscosity.

supercritical CO2； molecular dynamics； viscosity

2015-01-29

中国博士后科学基金资助项目，编号2014M552011；河南省基础与前沿项目，编号132300410104．

靳遵龙(1973-)，男，河南周口人，副教授，博士，主要从事过程装备与超临界流体强化传热研究，E-mail:zljin@zzu.edu.cn；通讯作者：王永庆(1979-)，男，河南林州人，讲师，博士，主要从事过程装备强化传热研究，E-mail:wangyq@zzu.edu.cn．

靳遵龙，李德雨，刘东来，等. 超临界CO2体系黏度分子动力学研究[J].郑州大学学报：理学版，2015，47(2)：101-103.

O552.3

A

1671-6841(2015)02-0101-03

10.3969/j.issn.1671-6841.2015.02.020