气体-表面相互作用中动量和能量分量间转化机制的分子动力学研究*

于航 张冉 杨帆 李桦

1) (国防科技大学空天科学学院, 长沙 410073)

2) (西北核技术研究所, 激光与物质相互作用国家重点实验室, 西安 710024)

1 引 言

随着微纳米技术和航天技术的高速发展, 稀薄气体流动相关的研究得到了越来越多的关注.对于处于稀薄环境中的微纳尺度系统和航天器来说, 流场中分子间碰撞频率的降低, 使得气体分子与固体表面之间的相互作用成为了影响流动状态的主要因素.然而由于其微观性和复杂性, 稀薄流动中气体-表面相互作用的机理并没有得到充分揭示[1,2],现有的模拟方法也不能准确地反映这一过程对流动的影响.因此无论从揭示物理现象本质的角度出发, 还是从为相关技术发展提供可信仿真数据的角度出发, 都有必要对稀薄流动中气体-表面的相互作用进行深入、细致的研究.

实际工程中采用的稀薄流动模拟方法大致可以分为两类, 一种是连续介质模拟方法, 其中包括附加了滑移边界条件的N-S 方程[3]、Burnett 方程[4]、Super-Burnett 方程[5]等, 这类方法在处理边界问题时, 都需要借助于滑移边界条件.另一类常用的方法被称为粒子类模拟方法, 如直接模拟蒙特卡罗法(DSMC)[6]和格子玻尔兹曼方法(LBM)[7], 为了描述气体分子与壁面之间的相互作用, 粒子类方法通常需要引入气体-表面相互作用模型, 即散射核模型.对于连续介质模拟方法来说, 现有的滑移边界条件, 无论是线性的还是非线性的, 都不能对壁面附近物理量的分布进行合理描述[8−11].很多研究人员也证明, 粒子类方法中广泛采用的Maxwell 和CLL (Cercignani-Lampis-Lord)散射核模型在模拟气体-表面相互作用方面也存在很多局限[12−14].例如, Woronowicz 和Rault[15]以及Bruno 等[12]对Xe 在GaSe 表面散射情况的研究表明, 这两种散射核模型都不能完全反映Xe的散射规律.并且,现有的散射核模型大多也不能描述气体分子在表面发生的吸附和捕获情况[16].从上述分析可知, 这两类实际工程中常用的稀薄流模拟方法都不能对气体-表面相互作用进行准确描述, 也不能充分反映气体-表面相互作用对流动状态产生的影响.为了修正和改进现有稀薄流动的模拟方法, 提高它们在模拟气体-表面相互作用方面的可信度, 有必要针对这一过程的特性和机理开展研究.

分子动力学(molecular dynamic, MD)[17]是一种从分子运动的角度出发, 直接根据分子间势能函数和牛顿第二定律对运动分子进行追踪的模拟方法, 可以准确地获得运动分子的全部参数, 在研究气体和表面的微观作用方面具有天然的优势.MD 方法被越来越多的学者应用到稀薄流动和气体-表面相互作用的有关研究中, 并且被证明可以获得与实验数据相符合的结果[18,19].Chirita 等[20]利用MD 方法研究了低速Ar 分子在Ni (001)表面散射过程中气体分子与固体表面之间的动量和能量交换, 并且获得了这一过程与分子入射角度、入射能量以及金属表面温度之间的依赖关系.Yamanishi 等[21]利用MD 方法研究了 O2, N2和Ar在石墨表面的散射, 并根据对模拟结果的分析, 提出了一种多级(multistage, MS)气体-表面相互作用模型.该模型被用到了DSMC 模拟中, 模拟结果可以很好地与分子束实验数据符合.曹炳阳[22]将MD 方法应用到了微尺度气体的模拟中, 并指出气体分子在固体表面附近的“俘获-逃逸”行为导致了氩分子在光滑铂表面上切向动量系数(TMAC)随温度的变化呈现指数衰减的关系.Finger 等[23]利用MD 方法研究了在气体-表面相互作用过程中, 气体分子的入射速度、晶体结构、分子间作用力和吸附层对气体分子散射和切向动量的影响.此外, 他们还根据气体分子与壁面的碰撞次数和反射速度方向对气体分子的散射方式进行了分类, 分析了不同散射类型对气体分子切向动量变化的影响.Ozhgibesov 和Leu 等[24−26]使用MD 方法研究了Ar 在光滑及粗糙W (110)表面的散射特性, 并基于模拟得到的反射速度和角度等信息构建了一种新的三维氩气-钨相互作用模型.Bao 等[27]利用MD 原理模拟了纳米管道内气态分子的速度滑移,结果表明, 在金属表面存在着一个气体分子的吸附层, 而这个吸附层的密度会随着温度的降低和气体-表面相互作用强度的增强而增大.他们还发现,管壁处的速度滑移会随着温度的升高和作用强度的减小而线性增加.宋诚谦[28]在一个利用GPU 加速的MD 模拟平台上计算了Ar 在Pt 表面的散射特性, 并利用得到的结果对CLL 模型进行了修正, 利用一个轴对称驻点流动的算例, 证明了修正后的CLL 模型可以实现与MD 方法较为一致的结果.

通过上述前人的工作可以发现, 稀薄条件下的气体-表面相互作用的影响因素众多, 并且对这一过程特性和机理的研究可以为相互作用模型的修正和改进提供有益的参考.动量和能量是反映气体分子运动状态的重要参数, 也是稀薄气体流动研究中的重要研究内容, 如Cao 等[29]曾对微/纳米机电系统中流-固界面的分子动量输运进行了充分而深入的总结, 综述了利用实验方法、理论方法和MD 模拟方法获得的分子动量输运的研究成果, 并详细介绍了流体黏性、极性、湿润度和壁面粗糙度等因素对这一过程的影响.前人关于分子散射后动量和能量的研究多集中于气体与壁面之间的动量、能量转化[20]和适应系数(MAC, EAC)[30,31], 却没有关注这一过程中气体分子动量、能量分量的变化.而根据我们之前的工作[32]可以发现, 在与壁面发生碰撞之后, 气体分子动量分量和能量分量之间的变化并不是一致的, 并且动量的两个分量对壁面粗糙度变化的敏感性也有着明显的差异.为了能够更深入地了解气体-表面的碰撞过程, 本文将利用三维非平衡MD 方法对Ar 分子与金属Pt 表面相互作用过程中动量和能量分量之间的变化和转化进行研究, 并考虑入射速度、角度和壁面粗糙对这一过程的影响.

2 物理模型及模拟细节

本文利用“Phantom”模型来模拟Pt 表面的原子结构[33], 这是一种基于Langevin 原理构造的恒温壁面模型.该模型中包含了三个原子层, 即利用虚拟弹簧力和势能函数模拟壁面原子与气体分子相互作用的真实原子层, 用来将第一层原子受到的作用力向底层传递的缓冲层, 以及稳定整个壁面系统的基底层.利用“Phantom”模型构建的金属表面是绝对光滑的, 而在工程应用中, 除了少数经过高精度加工的表面外, 绝大多数的表面都是粗糙的.因此, 为了研究壁面粗糙度对气体-表面相互作用过程中动量、能量分量之间的转化关系, 遵循文献[32]中的方法对“Phantom”表面模型进行修改, 使其更符合工程实际.本文选择了三种表面粗糙度构型, 即R = 0.0 Å, R = 0.5 Å和R = 1.0 Å, 结构如图1 所示.第一种粗糙度构型对应的是光滑表面,后两种粗糙度结构的原子堆高度分别为0.27 和0.54 nm.

图1 三种表面粗糙度构型Fig.1.Three surface roughness structures.

采用截断LJ-12-6 势能函数[34]来模拟气体分子与壁面之间的相互作用, 其具体形式如(1)式所示:

为了方便计算及分析, 将壁面原子的质量及势能参数设置为与气体分子一致, 即 σAr-Pt=σAr-Ar,εAr-Pt=εAr-Ar, 其中 σAr-Ar=3.405×10−10m 表示氩原子的直径, εAr-Ar=1.67×10−21J 为势能参数,rij代表分子间的中心距离, rc为截断半径.相比于原始LJ 势能函数, 截断半径 rc的引入可以在不影响计算精度的基础上减小计算量, 并且该形式的势能函数可以避免在截断半径处出现作用力不连续的情况.通过势能函数计算了气体分子与壁面的相互作用力之后, 为了兼顾稳定性和计算精度, 气体分子运动方程的积分通过velocity-Verlet[35]算法来实现.

在MD 方法中, 模拟单个分子入射有两种实现途径[36], 其中一种被称为“采样法”, 在这种方法中, 气体分子的速度满足特定温度下的Maxwell-Boltzmann 分布, 因此多被用来模拟自由来流条件下的气体-表面的相互作用[37].而为了能够模拟以特定速度和角度入射的气体分子, 本文采用的是分子束方法(molecular beam method)[38], 这种方法不仅降低了模拟系统的复杂度和不确定性, 还可以减小模拟过程的计算量.

本文采用的坐标系定义如图2 所示, 其中, 入射速度与Z 轴的夹角 θ 即为气体分子的入射极角,气体分子离开金属表面时的速度 vr可以分解为平行和垂直于XOY 平面的切向散射速度 Vt和法向反射速度 Vn.

图2 气体分子的入射极角 θ 和反射速度的分解在笛卡尔坐标系下的示意图Fig.2.Schematic diagram of the incident polar angle θ and the decompositions of molecular scattering velocity in Cartesian coordinate.

在进行MD 模拟时, 为了方便计算, 同时为了减小由于不同物理量的量级差异而产生的舍入误差, 需要对真实物理量进行约化处理.本文采用的约化基本单位是氩原子的LJ 势能参数, 包括直径参数 σ , 势阱 ε 和质量m.利用这些基本物理量, 在MD 模拟中的速度单位被重新定义为“ σ /τ ”, 其中代表约化时间.本文选取的分子的入射速度(Vi)都是特征温度下气体分子最可几速度其中 kB代表玻尔兹曼常数, 特征温度T取值为298 K)的倍数, 例如 Vi=4.0 即表示入射速度的数值是分子最可几速度值的4 倍.

3 结果与讨论

由于气体分子的动量是分子质量与速度的乘积, Ar 分子的质量为定值, 因此本章利用气体分子速度的值及变化代表其动量的大小与变化.我们在之前的工作中[32]详细研究了低速气体分子以45°的极角入射到金属表面上的动力学特征, 并且发现: 1) 在与壁面相互作用之后, 相同粗糙条件下,随着入射速度的增加, 切向和法向能量基本保持相等, 而两个动量分量之间的差距却会逐渐加大, 说明气体分子在金属壁面散射后动量和能量分量之间的变化和转化并不是一致的; 2) 切向动量和法向动量对壁面粗糙度的敏感性不同, 例如, 当入射速度为8.0 时, 气体分子在光滑表面上反射后会损失35%的切向动量, 法向动量的损失是16.5%,随着壁面粗糙度增加到1.0 Å, 切向动量的损失率达到了99.4%, 而法向动量的损失仅为25%.我们后续对高速入射时的散射特性的研究(结果汇总在表1 中, Vt和 Vn的定义同上文, Et和 En分别表 示散射后气体分子的切向能量和法向能量, E 表示散射分子的总能量.与文献[32]相同, “Re”表示反射表面的状态, 表中的散射数据是10 万个以相同条件入射的Ar 分子散射后速度和能量值的平均)也印证了上述现象.上述现象的出现说明, 气体分子在壁面发生碰撞的过程中, 动量、能量分量之间存在着复杂的转化机理, 有必要对其进行深入的研究.为了能够清晰地显示动量、能量在切向和法向之间的转化关系, 本章选取的气体分子的入射极角分别是法向分量占主导地位的 5°和切向分量占主导地位的 7 5°.

3.1 小极角入射时的动量、能量转化机制分析

首先对入射角为 5°时的算例进行分析, 气体分子在不同表面反射前后的速度和能量信息汇总在表1 中.可以看出, 当分子的入射方向与Z 轴夹角为 5°时, 气体分子速度和能量的切向分量相比于法向分量几乎可以忽略不计.而在气体分子与金属表面发生碰撞之后, 切向速度会进一步减小, 根据表面粗糙度的不同, 减小幅度为50%—94%不等.但是气体分子的切向能量却在碰撞之后增加了, 例如, 在入射速度为1.0 时, Ar 分子在光滑表面上反射后, 切向能量增加了近59 倍.由于气体分子接近和离开壁面的过程中, 壁面力场的做功会相互抵消, 所以这部分增加的能量可以看作是从法向转移来的.通过对散射后切向和法向能量值的分析发现, 当入射速度不小于2.0 时, 在某种特定的壁面条件下, 切向能量与法向能量的比值都只是在一个很小的范围内波动.例如, 对于光滑平面来说, 当入射速度在2.0—16.0 之间变化时, 切向能量与法向能量的比值维持在0.262—0.282 之间, 对于R =0.5 Å的平面, 这个范围是0.273—0.299, 当粗糙度进一步增加到1.0 Å时, 切向和法向能量之比稳定在 0.320—0.334 之间.同时, 还可以发现, 气体分子在固体表面散射后的总能量发生了损失.根据表1 中数据可以看出, 当气体分子以接近垂直的角度入射时, 粗糙度对动量、能量变化的影响不明显,因此, 为了减小计算量和简化仿真过程, 后续的相关计算都是在光滑表面条件下进行的.

表1 入射极角为 5° 时, 不同速度的气体分子在不同表面上散射后的平均速度和能量Table 1.When the incident polar angle is 5° , the mean velocity and energy of gas molecules that with different velocities after scattering on different surfaces.

为了对这一过程中动量、能量转换的机理进行更细致的研究, 进一步从微观角度对气体分子与表面碰撞过程中速度与能量的变化进行分析.图3是一个速度为2.0 的气体分子与表面碰撞过程中速度及能量随时间变化的曲线, 由于小角度入射时, 分子能量在壁面力作用区域的变化较为平稳,为了更清晰地反映能量的变化, 图3(b)只给出了气体分子与壁面碰撞前后一小段时间内的能量变化.需要注意的是, 由于表2 中数据是10 万个气体分子在表面散射后速度和能量值的平均, 而图3反映的是单个入射分子速度随时间的变化, 因此数值上会与表中数据有所差异.可以看出, 在以接近垂直的角度与金属表面碰撞的过程中, 气体分子法向速度的值会逐渐减小, 当减小为0 时, 意味着气体分子在这一方向的动能完全转化为了势能.由于受到壁面强大的排斥力, 在与壁面作用之后, 气体分子的法向动量会向相反的方向增加, 最终势能又转化为了分子的动能.在这一过程中, 法向的动量和动能会发生明显的损失, 而切向动能会在气体分子接近和离开壁面的过程中逐渐增大.

为了考察与壁面的相互作用过程中, 气体分子法向动量的损失及其与表面的适应程度, 进一步对分子束的法向动量适应系数(NMAC)进行研究,结果如图4 所示.NMAC 的定义如(2)式所示, 式中P 表示分子的平均动量, 下标“i”和“r”分别代表入射和反射,Tw表示壁温.

图3 入射角为 5° 时, 气体分子与光滑表面碰撞中速度和能量随时间的变化 (a)速度随时间的变化; (b)能量随时间的变化Fig.3.When the incident angle is 5° , the variations of molecular velocity and energy in the collision with a smooth surface: (a) Variation of molecular velocity; (b) variation of molecular energy.

图4 以 5° 极角入射时, 气体分子在不同粗糙表面的法向动量适应系数(NMAC)Fig.4.Normal momentum accommodation coefficient(NMAC) of gas molecules on different rough surfaces when the incident polar angle is 5°.

从图4 可以看出, 对于不同的表面,Vi=1.0时, 气体分子的NMAC 都保持在了一个较大的数值, 说明此时气体分子的法向动量基本上与表面达到了完全适应, 法向动量的损失明显.随着入射速度的进一步增加, 分子束的NMAC 发生了明显降低, 并会随着速度的增加而缓慢增大.同时还可以发现, 相同入射条件下, 金属表面的粗糙对气体分子的法向动量适应系数的数值影响不大.图5 为不同粗糙表面条件下, 动能的损失率随入射速度的变化, 整体来看, 动能的损失率会随入射速度的增加而增大.这是因为, 当气体分子的入射速度增加时,分子进入表面力场的深度就更深, 因此气体-表面之间的适应程度会更高, 这也就导致了气体分子的能量损失率的增大.而粗糙度的增加会使分子在金属表面的吸附概率增加[32], 尤其是分子的速度较低时吸附会变得明显.所以当 Vi=1.0 , R = 1.0 Å时, 气体分子的能量损失较大, 导致了此时的变化规律与整体趋势不一致.随着入射速度的增大, 气体分子在金属表面的吸附不再是影响其能量损失的主要因素.从图5 还可以看出, 在速度较高时,表面粗糙度对能量损失率的影响并不明显.

3.2 大极角入射时的动量、能量转化机制分析

图5 以 5° 极角入射时, 气体分子在不同粗糙表面的能量损失率Fig.5.Energy loss rate of gas molecules on different rough surfaces when the incident polar angle is 5°.

气体分子以 5°入射时, 动量和能量的法向分量占主导地位, 为了进一步研究切向动量和能量向法向的转移机制, 本节中将气体分子的入射极角设置为 7 5°.表2 列出了入射角为 7 5°时气体分子散射前后速度和能量的变化.可以看出, 与垂直入射时不同, 大极角入射时粗糙度对动量和能量转移的影响变得明显了.例如, 当速度为4.0 的粒子入射到光滑平面时, 反射后的切向和法向动量的绝对值都与入射前很接近, 只是法向动量的方向发生了改变,切向、法向能量的变化和总能量的损失也维持在2.5%—4.0%以内.而当表面的粗糙度为0.5 Å时,切向动量和能量的减小、法向动量和能量的增加以及总能量的损失就不能忽略了.当粗糙度继续增加到1.0 Å时, 切向动量的损失达到了90%以上, 法向能量几乎增加了4 倍, 总能量的损失也达到了34.5%.上述结果表明, 表面粗糙度会促进切向动量和能量向法向的转移.为了进一步探究气体分子各个动量和能量分量之间的转化关系, 分别对速度为2.0 的气体分子与光滑(R = 0 Å)和粗糙(R = 1.0 Å)表面碰撞过程中速度和能量随时间的变化进行研究, 结果分别如图6 和图7 所示(由于大角度入射时, 在整个壁面力场气体分子的能量都会有较为明显的波动, 只显示碰撞前后一小段能量变化的方式已不再适应, 因此本节给出的都是气体分子在整个壁面力作用区域速度和能量的变化).从图6 可以看出, 当气体分子以大极角入射到光滑壁面时, 散射规律非常接近于Maxwell[39]所假设的镜面反射, 具体表现为入射前后切向速度和能量的变化都很小, 并且这种可以维持入射前切向动量、能量的现象会随着入射速度的增加而更加明显.结合表2 中数据可以发现, 当 Vi=8.0 时, 入射气体分子的切向速度和切向动能分别为 1 7.293σ/τ 和149.53ε,在光滑表面反射后, 气体分子的切向速度和切向动能分别为 1 6.988σ/τ 和 1 46.57ε.同时, 气体分子的法向速度和法向能量在反射前后也变化不大, 只是法向速度的方向发生了改变, 这说明对于光滑表面, 入射速度的增大基本不会对动量、能量分量之间的转化产生影响, 气体分子的法向速度和能量变化只能依靠与壁面的热交换.与图6 中气体分子速度和能量的变化基本保持独立的现象不同,表面粗糙度的出现会改变气体分子与表面之间的动量、能量交换模式.如图7 所示, 当金属表面粗糙时, 在与表面碰撞之后, 分子的切向速度和能量会减小, 而法向速度和能量会增加, 说明粗糙度的存在促进了分子切向动量和能量向法向的转移.

表2 入射极角为 7 5° 时, 不同速度的气体分子在不同表面上散射后的平均速度和能量Table 2.When the incident polar angle is 7 5° , the mean velocity and energy of gas molecules that with different velocities after scattering on different surfaces.

图6 入射角为 7 5° 时, 气体分子与光滑表面碰撞中速度和能量随时间的变化 (a)速度随时间的变化; (b)能量随时间的变化Fig.6.When the incident angle is 7 5° , the variations of molecular velocity and energy in the collision with a smooth surface: (a) Variation of molecular velocity; (b) variation of molecular energy.

图7 入射角为 7 5° 时, 气体分子与粗糙表面碰撞中速度和能量随时间的变化 (a)速度随时间的变化; (b)能量随时间的变化Fig.7.When the incident angle is 7 5° , the variations of molecular velocity and energy in the collision with a rough surface: (a) Variation of molecular velocity; (b) variation of molecular energy.

为了研究不同入射条件下, 气体分子切向动量与表面的适应情况, 按照(3)式对大角度入射时的切向动量适应系数(TMAC)进行了计算, 结果如图8 所示, (3)式中 Pit和 Prt分别为气体分子的入射切向动量和反射切向动量的平均值.

由图8 中的数据可知, 对于光滑表面, 气体分子的TMAC 随着入射速度的增大而减小, 直至接近于0, 这说明当入射速度足够大时, 切向动量不会与表面进行明显的动量交换, 也就不会发生显著的动量转移了, 与上文得出的气体分子在光滑表面的散射基本满足镜面反射的结论相符.当表面粗糙度为0.5 Å时, 气体分子切向动量与表面的适应程度大幅提升, 但仍表现出随着入射速度增大而减小的趋势.随着表面粗糙度增加到1.0 Å, 不同入射速度下分子的TMAC 都接近于1, 此时气体分子的切向动量会与表面发生充分的交换, 这为切向动量向法向的转移创造了条件.图9 为不同表面条件下, 能量的损失率随入射速度的变化.可以看出,动能损失率都会随入射速度的增加而增大, 但当速度增大到一定值后, 能量损失率就不再继续变化了, 与入射极角为 5°时不同, 入射角为 7 5°时, 粗糙度的值会对气体分子在表面上的能量耗散产生显著的影响, 粗糙度越大, 总动能的损失也就越大.

图9 以 7 5° 极角入射时, 气体分子在不同粗糙表面的能量损失率Fig.9.Energy loss rate of gas molecules on different rough surfaces when the incident polar angle is 7 5°.

4 结 论

本文对两种入射角度下气体分子动量和能量分量之间的转化进行了研究, 并考察了壁面粗糙度和入射速度对转化机制的影响.当以接近垂直的角度与壁面碰撞时, 气体分子切向和法向的动量都会发生损失, 而法向的能量会向切向转移, 并且当分子速度较大时(≥ 2.0)其转移的比率受分子入射能量和表面粗糙度的影响不明显.气体分子在散射后总能量的损失会随入射速度的增加而增大, 而不会随着表面粗糙度的变化产生明显差异.当气体分子以 7 5°入射时, 壁面的粗糙度会对分子动量和能量的转化机制产生显著的影响.在与光滑壁面碰撞之后, 气体分子的动量和能量值都基本保持不变, 只是动量的方向发生了反转, 其运动状态接近于镜面反射, 动量、能量分量之间的转化不明显.而粗糙度的引入则使气体分子与表面之间的适应程度增强, 促进了分子切向动量和动能向法向的转移.大角度入射时, 分子总能量的损失对入射速度的变化不敏感, 但是会随着粗糙度的增大而发生明显的增加.本文所做的工作不仅是对气体-表面相互作用机理的部分揭示, 还可以为稀薄气体流动和气体-表面相互作用的高保真模拟提供有益的参考.