量子流体氦黏度的分子动力学模拟

　　引　言

　　黏度是流体的基本迁移物性[1],通常用动力黏性系数η来表征单位速度梯度作用下的切应力,它反映了黏滞性的动力性质[2]。近年来,计算机技术的发展促使部分有关流体的问题可以通过计算机模拟手段予以解决。鉴于传输参数黏度的重要性,很多学者都针对氩这种典型的LJ (Lennand-Jones)流体进行了关于黏度的MD模拟研究[3-5]。但在处理一些较轻原子或分子的平动或转动,或振动频率ν满足hν>kBT的振动时,需要考虑量子效应[6]。氦就是量子流体的典型代表,且随温度的降低,其量子效应增强。传统LJ势函数并未考虑到量子效应,采用此模型模拟氦的黏度,势必产生一定的误差。因此,在LJ模型的基础上,考虑对其进行量子效应修正是十分必要的。

　　1　势能模型

　　分子势能模型的选择会直接影响模拟的结果[7]。在带有量子效应修正的分子动力学研究领域,还较少涉及传输参数的模拟[8]。目前可以通过FH (Feynman-Hibbs)势函数和WK (Wigner-Kirkwood)势函数实现LJ流体的量子效应修

　　正[9]。保留二次波动项的势函数被称为QFH(quadratic Feynman-Hibbs)函数,将LJ函数表达式代入QFH函数,并将其写成量纲1的形式,得到

　　式中　β=1/ TkB, J-1,kB为Boltzmann常数; =h/2π,h=6·6256×10-34J·s;μ=m/2, kg。从QFH势能函数的表达式可知,其量子校正项是温度的函数,温度越低量子效应越强,即对LJ函数的修正性增强,如图1所示。

　　有些学者利用带有量子效应的势函数实现了对量子流体氦的MD研究,获得了流体黏度等其他传输参数,但所得到的状态点非常少[4]。粒子之间的作用力可通过式(2)计算[10]

　　由此可得到量纲1形式的QFH势函数作用力计算方程式(3)和式(4)

　　式(3)和式(4)的第二项是量子校正项,系数a随温度的变化关系如图2所示。

　　2　模拟原理与细节

　　固体壁面和流体均由简单的球形分子组成,分子运动遵循牛顿第二运动定律。对于流体,分子i的速度.ri和动量pi由式(6)给出[11]

　　式中　m是分子的质量,fij表示发生在分子i与j之间的作用力,而Fe是指由于外界力场的作用产生的力。固体壁面通常用两层壁面分子表示,质量与流体分子的质量相同,运动方程与式(6)相似,只不过将弹簧的弹力替代了式(6)中的外力[11]。模拟中,采用外加重力力场的方法驱动流体产生流动。作用在分子上的重力力场使其在X方向产生流动,而能调整的就是在模拟中所能采用的重力加速度常数。外界力场驱动流体流动,同时也对流体做功。部分功转化为热量,使流体和壁面的温度升高。为保证系统温度始终处于某一固定温度,需要启动一些具有特殊功能的子程序,本模拟中采用的方法是在模拟进行的每一个时间步,都对壁面温度进行调整,因此壁面可以维持在一个固定的温度。流体和壁面分子的初始构型由预处理程序事先得到,并通过设置输入文件将其数据输入到模拟程序中。有些学者在处理流体黏度时将模拟得到的速度分布根据带有无滑移边界条件的Navier-Stokes的等温解进行拟合,从而可以利用牛顿流体力学的理论和关系式分析讨论模拟的结果[11-12]。根据经典流体力学理论中关于黏度η的物理理解:当速度梯度du/ u/dx/ x=1时( u和x分别表示u和x的单位),τ=η,即η表征单位速度梯度作用下的剪切应力。因此,可以写出剪切应力与速度分布的关系式