双原子氧气液界面特性的分子动力学模拟

毛志红 包福兵 余 霞

(中国计量学院计量测试工程学院 杭州 310018)

1 引言

气液界面热力学行为一直是相变传热传质研究的重点。由于气液界面层很薄，并且界面不固定，使得采用理论分析和实验方法研究界面特性十分困难。近些年来，越来越多的学者采用数值模拟，特别是采用分子动力学方法来研究气液相变界面特性［1］。

Trokhymchuk 和 Alejandre［2］采用分子动力学方法计算得到了Lennard－Jones流体氩气液相平衡时的密度、法向与切向压力张量分量以及界面张力等的分布曲线，总结了气液相平衡时气液相和界面的密度、压力张量分量等的变化规律。Frezzotti［3］等人采用分子动力学方法系统建立了气液界面模型，指出非平衡态时，在气体与气液界面接触处形成努森层，描述了纯质和多组分流体蒸发凝结时努森层的作用，并分析了边界条件对气液界面性质的影响。Lee［4］采用分子动力学方法，从基于维里定理、Irving－Kirkwood模型和Harasima模型的计算方法出发，模拟得到了氩流体在气液平界面处的压力张量，并对不同方法得到的结果进行了对比，发现采用这3种方法模拟得到的切向压力张量分量和表面张力基本相同，但通过Irving－Kirkwood模型得到的法向压力张量分量在低温时几乎不变，与另外两者有很大不同。最近，Chilukoti［5］等人对丁烷、己烷、癸烷等在界面附近的自扩散和界面结构进行了分子动力学模拟，发现在温度相同时，这3种烷在界面处的排序和分子链形类似。

中国国内学者对气液界面微观性质的研究开始得较迟。过增元院士研究小组［6］首次采用分子动力学模拟方法研究了相界面的形成和变化问题，提出气液界面过渡区的宽度与系统的温度有关，并研究了简单氩流体的气液界面的微观特性，分析了界面区的热力学性质，得到气液界面的密度、压力张量分量及局部温度的分布等。陶文铨院士研究小组［7］研究了气液界面表面张力的大小与系统尺度的关系，分析了界面温度和入射分子数对界面凝结过程的影响。曾丹苓教授研究小组［8］系统研究了纳米流体的输运特性，并探究了纳米液滴在蒸发过程中界面现象的物理机理。

但是到目前为止，中国国内外关于气液界面性质的研究大多针对单原子分子，而实际工程中的气体，往往是多原子分子，例如空分精馏过程中氧气和氮气，对这类分子的研究较少。即使有，也多是采用单原子模型来研究，这种模拟方法与实际分子结构存在一定误差。氧气作为典型的双原子分子，具有代表性。随着钢铁、化工、国防、医疗的发展，国内各行业对氧气的需求量越来越大。氧气主要采用深冷法制取，为提高工业生产效率，改善生活，对氧的低温热力学性质的研究提出了迫切需求。本文采用分子动力学方法，从Lennard－Jones势能函数模型出发，考虑氧分子中原子间键的作用，研究双原子氧的气液界面性质。

2 模拟方法

采用分子动力学方法来模拟双原子氧的气液界面性质，氧分子由2个氧原子通过键连接构成。氧原子间的势能采用平移后的Lennard－Jones 12－6势能函数模型，其表达式为:

式中:ε为势阱深度;σ为原子作用直径;rij为i和j两原子间距;rc为原子作用截断半径，在截断半径处，两原子间势能为零。

氧分子在研究中一般被看作是一种具有弹性的分子，氧分子中原子间键的作用常常采用一种谐振势函数表示:

式中:K为弹性常数;xe为平衡键长，氧分子内键长［9］xe=1.208 Å。

基于NVT系综，计算区域如图1所示，左右两侧为气相区，中间区域为液相区，系统尺寸为55Å×165Å×55Å，3个方向都采用周期性边界条件。为了提高计算效率，采用元胞链表法来计算粒子间作用力，采用Verlet蛙跳法积分粒子的位置和速度。计算中，粒子数N=3 000，时间步长取为Δt=1 fs。一开始所有分子均匀分布在液相区域，根据温度设定初始速度。根据文献结果及对计算数据的分析，系统经过20万步时达到平衡，再统计之后的180万步的结果进行分析。

为了统计界面性质，沿y方向把计算区域平均分成Ns=165层。根据每层中原子数和速度，统计得到每层体相的密度、饱和压强及表面张力，获得平衡态气相饱和压强、气相及液相密度、界面区表面张力等性质。每层厚度为Δy的分子数密度统计公式:

图1 气液相变系统示意图Fig.1 Schematic diagram of gas-liquid phase change system

式中:Vs=Lx×Δy×Lz。

根据Kirkwood－Buff对表面张力的定义，结合本文建立的气液界面模拟体系，表面张力γ可以通过式(4)来计算:

其中:PN和PT分别为界面处的法向压力张量分量和切向压力张量分量。

式中:［］表示求总体均值，ρ(y)、T(y)、V(y)分别表示局部数密度、温度和体积;dφ(r)/d r为势函数对r的导数;kB为波尔兹曼常数，1.38×10－23J/K。

3 氧原子参数确定

采用分子动力学方法模拟气液界面性质过程中，势能参数对计算结果影响很大。以往有关气液界面性质研究中，通往把氧分子看成一个整体，忽略了分子中原子的振动和取向对整体的影响，而采用双原子模型研究氧气液界面性质的文献并不多见。Arora［9］采用双原子模型研究了氧分子在碳纳米管中的运动中，使用了如下的参数:ε/kB=52 K、σ =2.99 Å、xe=1.208 Å、K=11.77 N/cm，该参数能够较好地描述氧分子与碳纳米管的相互作用。但是，当采用上述双原子模型和参数模拟存在相变现象的氧气的气液界面时，所得结果与实验数据［10］相比误差较大。

分子间的引力和斥力作用对气液相变影响巨大。引力越大，界面中的分子越难摆脱液相区分子的束缚进入气相区，而气相区分子越容易被界面中的分子捕捉而进入液相;斥力越大，界面中的分子越容易受到液相区分子的斥力作用而成为气相分子，而气相区分子越难进入界面区进而进入液相。势阱深度ε代表了分子间作用力的强弱，对气液界面性质有较大影响，因此，为简化计算，模拟中，根据实验测得的气液相密度，调整势阱深度ε，而保持其他参数不变。根据典型温度T=80 K，截断半径R=4.5σ，对势阱深度 ε/kB=52、57.97、59.36、59.94、60.17、60.29、61.6 K共7种情况下模拟结果与实验数据(ρg=1.468 4 kg/m3，ρl=1 190.5 kg/m3)［10］进行了比较，具体见表1。

表1 80 K时，氧在不同ε下气液相密度与实验值的比较Table 1 Gas-liquid density of oxygen compared with experimental value under differentεat 80 K

图2 氧气气液相密度随温度的变化Fig.2 Gas-liquid density of oxygen under different temperatures

图3 氧气的饱和压强随温度的变化Fig.3 Saturation pressure of oxygen under different temperatures

从表1中可以看出，随着ε/kB的不断增大，ρg不断减小，ρl不断增大。当 ε/kB=60.17 K时，气液相密度误差在2.5%以内，与实验值十分接近，因此，在后续计算中都采用这个值。图2和图3分别给出了系统达到平衡时氧的气液相密度随温度的变化以及氧的饱和压强随温度的变化曲线，并与实验值进行了比较。从图中可以看出，随着温度升高，液相区密度逐渐减小，气相区密度逐渐增大。计算中，当温度达到154.58 K时，系统密度处于均匀状态，定义此时温度为氧的临界温度，模拟值与实验值符合很好，验证了双原子氧模拟程序的正确性以及确定的势能参数的可靠性。

4 结果及分析

根据确定的氧原子势能参数，采用分子动力学方法模拟获得不同温度下氧的压力张量分量分布曲线，图4给出了60、90、120和140 K时的气液界面的法向压力张量分量PN和切向压力张量分量PT在界面附近的分布。从图中可以看出，法向压力张量分量和切向压力张量分量在气相区域和液相区域的值基本相等，液相区域压力张量分量的波动比气相区域的波动较大，而在气液界面处存在较大的势垒和势阱，与此相比较，气液区域波动可忽略不计。法向压力张量分量PN先下降再上升，最后在平衡位置波动。切向压力张量分量PT先下降再在平衡位置波动。

图4 不同温度下氧气压力张量分量分布曲线Fig.4 Distribution curve of oxygen stress tensor component under different temperature

根据式(4)对表面张力的定义，可以得到界面附 近的表面张力分布及界面区的表面张力。图5给出了不同温度下界面附近氧气表面张力的分布。由图可以看出，在气相和液相区域，不存在界面，此时表面张力的值为0。而在气液界面区，表面张力值较大，表面张力的值随着温度的升高而减小。图6给出了不同温度下氧气的表面张力值，并与实验值［11］进行了比较。从图中可以看出，模拟值与实验值符合的很好。表面张力随温度升高而线性减小。当温度达到临界温度154.58 K时，表面张力值减小到0，不存在界面区。

图5 表面张力在界面附近的分布Fig.5 Distribution of surface tension near the interface

图6 表面张力随温度的变化Fig.6 Variation of surface tensions with temperatures

4 结论

采用Lennard－Jones势能模型，对双原子氧的气液界面特性进行了分子动力学模拟研究。通过计算结果与实验值的对比，发现当氧原子的ε/kB=60.17 K时，计算得到的氧气的气液相密度以及饱和压强随饱和温度的变化曲线都与实验结果符合得很好。通过计算发现，气相和液相区域的切向和法相压力张量分量相等，这些区域中的表面张力为零。界面区压力张量的法向与切向分量存在较大的势阱和势垒，界面区的表面张力随温度线性降低，当达到临界温度时，表面张力减小到0。通过模拟值与实验值的对比，验证了双原子模型及氧原子势能参数的正确性，为进一步研究氧原子在各种氧化物中的性质提供了帮助。

1 陈进，徐征，王惠龄.Ar－Kr低温界面传热过程分子动力学模拟［J］.低温工程，2005(3):45－46.

2 Trokhymchuk A，Alejandre J.Computer simulations of liquid/vapor interface in Lennard－Jones fluids:Some questions and answers［J］.Journal of Computational Physics，1999，111(18):8510－8523.

3 Frezzotti A.Boundary conditions at the vapor－liquid interface ［J］.Chem.Phys.，2011，23(3):030609－1－030609－9.

4 Lee SH.Pressure Analyses at the Planar Surface of Liquid－Vapor Argon by a Test－Area Molecular Dynamics Simulation［J］.Bulletin of the Korean Chemical Society，2012，33(9):3039－3042.

5 Chilukoti H K，Kikugawa G，Ohara T.A molecular dynamics study on transport properties and structure at the liquid－vapor interfaces of alkanes［J］.International Journal of Heat and Mass Transfer，2013，59:144－154.

6 王遵敬，陈 民，过增元.汽液界面动力学行为与热力学性质的分子动力学研究［J］.工程热物理学报，2001，22(1):70－73.

7 熊建银，陶文铨，何雅玲.气液界面特性和凝结过程的分子动力学研究［C］.中国工程热物理学会第十一届年会，北京:2005.

8 蔡治勇，刘娟芳，刘 朝，等.纳米液滴蒸发过程的分子动力学模拟［J］.工程热物理学报，2008，29(10):1630－1632.

9 Arora G，Sandler SI.Air separation by single wall carbon nanotubes:Mass transport and kinetic selectivity［J］.Journal of Chemical Physics，2006，124:084702.

10 陈国邦，黄永华，包 锐.低温流体热物性［M］.北京:国防工业出版社，2006.

11 Lemmon E W，Penoncello S G.The Surface Tension of Air and Air Component Mixtures［J］.Advances in Cryogenic Engineering，1994，39:1927－1934.