张烨 张冉 常青 李桦

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

1 引 言

近年来,随着微纳机电系统(MEMS/NEMS)的快速发展,纳米尺度气体流动问题受到广泛关注,逐渐成为一个重要的研究领域[1-4].纳米尺度气体流动具有广阔的应用前景,比如生物样本提取、集成电路冷却以及流动主动控制[5].然而,纳米尺度下的流动特性与宏观流动相比有很大差异.为了理解纳米尺度流动特性,有必要对其机理展开深入研究.

气体的稀薄程度一般用克努森数(Kn=λ/L,其中λ为气体分子平均自由程,L为流动特征长度)来表示.微纳尺度气体流动中,分子平均自由程与流动特征尺度相比不再是小量,连续介质假设不再成立,气体表现出明显的稀薄效应.同时,尺度的减小会导致面积与体积比增大,与壁面相关的物理量对流动的影响更加突出.当Kn≤0.01 时,流动处于滑移流区域,Navier-Stokes (NS)方程在主流区域仍然有效,仅在壁面处出现速度滑移,附加具体的滑移边界条件仍然可以求解.Maxwell滑移边界条件[6]以及一些能够反映固体界面处流体流动特性的滑移模型[7-9]已经在工程中得到广泛应用,但是NS方程的线性本构关系决定了滑移边界条件对于壁面附近物理量非线性分布的描述不足.随着克努森数的增大,平均自由程和特征长度处于同一量级,流动处于过渡流区域,克努森层变厚,壁面效应影响更加突出,流场表现出强烈的非平衡和非线性特征,基于连续性假设的NS方程彻底失效.而描述稀薄气体的最完备的玻尔兹曼方程由于自变量维数高,数学形式复杂,求解起来十分困难.对玻尔兹曼方程高阶展开得到的近似方程比如Burnett方程、Super-Burnett方程,虽然较NS方程有一定优势,但仍然不能完全描述过渡区稀薄气体的强非线性关系和非平衡效应.且由于稳定性等诸多问题没有解决,表达形式复杂,物理意义不明确,高阶方程难以推广应用.因此,人们考虑采用数值模拟方法研究稀薄气体流动.直接模拟蒙特卡罗(direct simulation Monte Carlo,DSMC)方法[10,11]是以分子混沌假设为前提的一种基于分子模型的方法,能够有效模拟过渡流区域稀薄气体流动.但是DSMC方法在处理边界处气体与壁面的相互作用时,仍然需要提供动量和能量适应系数.对于不同的气体和壁面,其适应系数相差很大,因此该方法不适合研究壁面占主导因素的纳米尺度流动.

分子动力学(molecular dynamics,MD)方法[12]作为一种确定论的模拟方法,从原子运动的角度出发,根据牛顿第二定律和分子间势能函数追踪每个分子的运动,不仅能获得流场的流动细节,还能获得气体分子在壁面处的散射规律,因此在研究微纳尺度流动方面具有独到的优势.基于MD模拟,许多学者针对微纳通道流动开展了广泛的研究.曹炳阳等[13,14]采用MD模拟方法对滑移区微纳通道内的气体流动进行了研究,探讨了速度滑移与势能作用强度、系统温度以及壁面粗糙度之间的关系.Priezjev[15]采用MD方法研究了具有纳米尺度纹理的壁面构成的通道内的液体流动规律,发现近壁面第一层液体的有效滑移长度和耗散系数依赖于湿壁区域的面积.Sun和Li[16]采用非平衡分子动力学方法计算了光滑微纳通道中的切向动量适应系数和法向动量适应系数,通过等温伯肃叶流动分析了温度、气体-壁面作用长度以及气体分子的吸附现象对动量适应系数的影响.Spijker等[17]采用MD方法模拟了纳米通道内的气体流动,统计了气体分子与壁面碰撞前后的速度分布,并在分析入射与反射速度间的关系时引入了速度相关性云图.Xie和Liu[18]模拟了纳米通道中混合气体的Poiseuille流动,发现壁面作用力使得通道内气体混合物呈现出化学成分与物理结构不再均匀的特殊现象.Kamali和Kharazmi[19]用MD方法研究了不同的壁面亲疏水性和粗糙单元长宽比条件下,液氩纳米尺度伯肃叶流动的密度和速度分布.Barisik和Beskok[20,21]引入了一种智能壁面分子动力学方法,假设壁面原子静止,发现通道内近壁区域气体的密度、速度与压力的变化仅由壁面力场确定,不受通道高度影响的规律.Noorian等[22]应用MD方法研究了液氩在微纳通道中的流动,结果显示壁面吸引能和粗糙高度的增加,会加强壁面附近处的密度分层现象.Bao等[23,24]采用MD方法对有限长度纳米通道内的压力驱动流动进行研究,得到了压力差与温度和克努森数无关而与通道纵横比及壁面性质相关的结论.To等[25]用MD模拟方法研究了微纳通道中非平衡气体的平均自由程分布情况,发现平均自由程在壁面附近显著减小,而对流动类型不敏感.Liakopoulos等[26]采用非平衡分子动力学方法研究了壁面的亲疏水性和粗糙度对纳米通道内摩擦系数的影响.Lim等[27]用MD方法模拟了He,Ne,Ar,Kr和Xe气体的流动,研究了分子质量和吸附效应对切向动量适应系数的影响,发现切向动量适应系数随分子质量的变化并不单调.王胜等[28]采用MD方法研究了液体在非对称性浸润性粗糙纳米通道内的流动与传热过程,发现随着两侧壁面浸润性不对称程度的增加,流体反转温度分布更加明显.

从前人的研究来看,壁面的影响在微纳通道流动中占据着主导作用.为了进一步研究壁面效应对纳米尺度气体流动的影响,我们在之前的工作中研究了通道高度和壁面势能系数对纳米通道内气体剪切流动的影响[29],但是对于具有一定刚度和粗糙度构型的壁面,其对流场物理量的分布以及近壁区域流动特性的影响研究还不是很充分.实际上通道壁面的原子并不是静止的,而是存在着分子热运动.纳米通道的壁面通常用Einstein固体理论[30]来构造,该模型通过弹性系数来模拟壁面原子的热运动,尽管有一些学者研究了弹性系数对纳米通道内液体流动的影响规律[31-33],但是在大部分微通道气体流动研究中,由于壁面原子量巨大,出于节省计算量的目的,都将壁面原子设置为静止不动[20,21],这种处理无疑忽略了壁面原子振动对气体流动的影响,因此有必要研究壁面刚度对气体流动的影响.此外,除了少数高精度加工的壁面外,绝大多数的壁面都是粗糙的.对于纳米通道内的气体流动,粗糙度的尺寸相比于流动的特征长度并不是小量,粗糙度对气体分子与壁面的能量和动量交换有着更加突出的影响,且粗糙度的存在增大了壁面力场的穿透深度,因此研究粗糙度对近壁区域气体流动的影响具有重要意义.系统温度的变化对纳米通道稀薄气体流动的影响主要通过两方面来体现.一方面,气体温度的变化会影响通道内气体分子之间的碰撞频率,而碰撞频率的变化又会导致气体的黏性系数发生改变,进而对流动产生影响;另一方面,温度的变化会影响气体分子与壁面的动量和能量交换,从而导致气体分子对壁面的切向动量适应系数的变化.因此温度也是纳米通道流动中的一个重要变量.

本文采用三维MD模拟方法,针对处于过渡区的纳米通道内气体剪切流动展开模拟,改变壁面刚度、壁面粗糙度以及系统温度,获得通道内流场的密度、速度、正应力及剪切应力分布,以分析壁面条件和系统温度对流动的影响规律,并探讨近壁区域内物理量的变化,进一步计算和分析切向动量适应系数(TMAC)的大小和变化规律.

2 三维分子动力学模拟

2.1 物理模型及模拟细节

本文采用三维非平衡分子动力学方法来模拟微纳通道中气态氩Ar的剪切流动,其模型示意图如图1所示.单原子气体Ar分布在由上下两块无限长平板构成的高度为H=10.9 nm的纳米通道内,X,Y方向采用周期性边界条件,模拟区域的长宽尺寸为LX×LY=65.38 nm×65.38 nm.上下壁面以的速度沿相反方向运动,驱动气体做剪切流动,其中M为马赫数,本文模拟时取为0.2;γ为比热比(对于单原子气体,γ=5/3);玻尔兹曼常数kB=1.3806×10-23J·K—1;T为气体温度(默认情况下T取298 K);m为气体分子质量.对于单原子气体 Ar,其分子质量mAr=6.63×10-26kg,分子直径σAr=0.3405 nm,势能参数εAr=1.67×10-27J.为了方便计算与分析,本文中壁面原子参数与气体分子保持一致,即mwall=mAr,σwall=σAr,εwall=εAr.

图1 纳米通道内气体剪切流动示意图Fig.1.Schematic of the shear-driven gas flow in the nanochannel.

气体-气体和气体-壁面之间的相互作用采用截断Lennard-Jones (L-J)6-12势能函数[34]模拟,定义为

初始时刻,气体分子均匀分布在通道内,且按照Maxwell分布产生初始随机速度.速度积分采用Velocity-Verlet算法,时间步长Δ=0.002τ,系统运行 2×106步后达到稳定,稳定后继续运行较长时间步,以确保所有宏观量都达到热力学稳定状态.为了捕捉近壁区域的物理量变化,将模拟区域沿Z方向细分为若干等间距(0.2σAr)的统计层,以统计流场中各宏观量的细致分布.

2.2 壁面刚度模型

壁面原子按照面心立方晶格点阵(FCC)排列,晶格常数为 1.6σAr,采用Einstein壁面模型[35],使壁面原子在其初始晶格点位置以相同的频率做简谐振动,简谐振动由连接分子和晶格点的虚拟弹性力维持.其振动的弹性系数K表征了壁面的刚度,由Lindeman准则[36]〈取值〉,通常要避免壁面〈整体〉移动较大,限制条件为 dA2/d2＜0.023,其中dA2代表分子相对平衡位置振动时的均方位移,d为壁面中紧邻的两个分子之间的距离.根据壁面原子类型的不同,可以根据其对应的Einstein温度来确定虚拟弹簧力的弹性系数：

其中m为壁面原子的质量,TE为固体的Einstein温度,h为普朗克常数.通过这种方法确定的弹性系数既满足Lindemann准则,还可以模拟壁面原子热运动及原子结构,能够反映真实的壁面状况.本文在默认情况下取K=1000ε/σ2,壁面刚度影响问题研究的计算条件如表1所列.

2.3 壁面粗糙度模型

本文构造的粗糙度为一种典型的金字塔形状的粗糙度模型[37],其几何形状示意图见图2.由于构成粗糙壁面的原子均位于晶格点的位置上,因此该形状的粗糙度具有合乎物理实际的倾斜度,在结构上具有很好的稳定性.本文设计了三种不同大小的粗糙度单元,高度分别为0.27,0.54及0.81 nm.

图2 壁面纳米粗糙度构型Fig.2.Schematic diagram of the roughness geometry.

相比于光滑壁面,粗糙壁面的粗糙度具有不可忽视的物理尺寸,通道的特征长度会因此而发生一定的变化.图3为粗糙通道的有效高度和光滑壁面的高度计算示意图,在模拟过程中,为了保证各算例的克努森数保持相同,通道的高度必须一致.对于不同的粗糙壁面,统一地将粗糙通道的上边界和下边界分别放置在下壁面和上壁面粗糙度的中心位置.由于气体分子会进入粗糙度内部的间隙中,而在统计流场中的物理信息分布时,这部分气体分子也是需要考虑的,因此,就会导致在后续的图表中出现粗糙通道的真实高度(Z)为负值及约化高度(Z/H)小于—0.5的情况.粗糙度影响问题研究的计算条件如表1所列.

表1 计算条件设置Table 1.Parameters of the study on the influences of surface stiffness,roughness and temperature.

2.4 温度控制方法

所有的系统都会与外界发生相互作用,比如受到外界温度和压力的影响.进行MD模拟时通常需要构造一个等温或等压的系综,一般来说,模拟系统的分子数目和体积比较容易保持恒定,而温度和压力则不同.由于系统中的粒子不停地运动与碰撞,并伴随着动量、能量的交换,因此,温度和压力是不断变化的.为了使系统温度或压力保持恒定,有必要对温度及压力进行适当的控制和调节,也称作热浴.本文采用Nose-Hoover热浴方法[38],则模拟系统的运动方程可写作：

式中F为粒子i所受到的作用力;η为热浴阻尼系数,

合作式治理路径主张社区既非国家垄断公共权力的领域，也非完全由社会自治的场所，而是国家权力和社会权力共同作用、合作治理的领域。研究者认为，社区治理是最接近于“治理”原义的治理形态，它是由多元主体共同参与、管理、决策和协商的合作治理结构，具有基层草根性、利益相关性与主体多元性的特征。⑤中国的社区治理要完成从威权式治理向参与式治理体制转型，赋予居民参与社区公共事务的权力，让居民自己决定关系切身利益的社区公共事务，形成以社区需求为中心的权力和资源分配体制，发展社区居民自我治理的理念和能力，让社区发展主体——居民能够有效地介入到社区建设的决策和管理过程中。⑥

其中T0为目标 控 制 温度,T(t)为当前时刻的 温 度,通过阻尼系数使流场温度保持在目标温度下.本文在默认情况下温度取为298 K,温度影响问题研究的计算条件如表1所列.

3 结果与讨论

3.1 壁面弹性系数对流动的影响分析

为了研究壁面弹性系数对纳米通道内气体流动的影响规律,本节设置了一组算例,弹簧系数设置在 3 00ε/σ2—1 000ε/σ2之间,作为对比,同时计算了壁面原子固定的情况,也即弹性系数无穷大.

图4为纳米通道内的速度分布曲线,由图可知,壁面弹性系数对主流区的速度分布影响不大,但对近壁面区域的速度分布有一定影响.随着壁面弹性系数的增大,气体分子在近壁面区域的速度峰值 分 别 为 ：—51.17,—50.08,—49.72,—49.23,及—47.26 m/s.从原子运动的角度出发,当气体分子与壁面相互作用时,壁面原子的振动频率越大,则气体分子就越不容易在壁面停留,从而导致气体分子与壁面的动量适应性变差,气体分子在近壁面区域的速度峰值随之减小.

图4 不同壁面弹性系数条件下流场速度分布 (a)无量纲速度分布;(b)近壁面区域速度分布Fig.4.Velocity profiles with different surface stiffness：(a)Normalized velocity profiles;(b)near-wall velocity profiles.

图5—图7分别为不同壁面弹性系数条件下,纳米通道中的密度、正应力及剪切应力分布.由图可知,壁面弹性系数对流场内上述物理量的分布影响不是很明显.不同壁面弹性系数条件下,流场中心区域的密度和剪切应力如表2所列,通过剪切应力与克努森数的关系[39,40],得到不同壁面弹性系数下的TMAC分别为0.740,0.736,0.731,0.731,0.726,总体上呈减小趋势,但影响较小.

3.2 壁面粗糙度对流动的影响分析

为了研究粗糙度对纳米尺度气体流动的影响规律,改变通道壁面的粗糙高度,模拟了不同粗糙度大小情况下的纳米尺度气体流动,统计获得了流场内速度、密度、正应力及剪切应力分布.

图5 不同壁面弹性系数条件下纳米通道中的密度分布Fig.5.Density distribution with different surface stiffness.

图6 不同壁面弹性系数条件下纳米通道中的正应力分布Fig.6.Normal stress distribution with different surface stiffness.

图7 不同壁面弹性系数条件下纳米通道中的剪切应力分布Fig.7.Shear stress distribution with different surface stiffness.

图8为不同粗糙通道内气体的无量纲速度分布,可知无论是主流区域还是近壁面区域,粗糙度对流动的影响都非常明显.而且由于粗糙度的存在,在壁面出现了一个粗糙度影响区域,在这个区域内,也存在着一定的速度分布规律.为了更加细致地了解这一区域的速度分布情况,截取了距下壁面边界2.0 nm以下范围内的无量纲速度分布曲线,如图9所示.从图9可以看出,在主流区域,气体的速度随着壁面粗糙度的增大而增大,当粗糙度高度为0.54和0.81 nm时,主流区域的速度分布趋于一致.此外,由于壁面粗糙度的增大,主流区的线性分布特性也受到了一定的影响,线性区域的范围减小.在近壁面区域,随着粗糙度的增大,流体速度迅速增大并接近于壁面速度的大小,当粗糙度足够大(0.81 nm)时,几乎整个粗糙单元缝隙中的气体流速都与壁面速度相等.

表2 不同壁面弹性系数条件下纳米通道中的剪切应力及TMAC大小Table 2.Shear stress and TMAC values with different surface stiffness.

图8 不同壁面粗糙度下的流场无量纲速度分布Fig.8.Normalized velocity profiles with different surface roughness.

图9 不同壁面粗糙度下近壁区域流场无量纲速度分布Fig.9.Near-wall normalized velocity profiles with different surface roughness.

图10为不同粗糙度高度下近壁面区域的密度分布曲线,可知随着壁面粗糙度的增大,气体的在近壁面附近的峰值随之增大,同时粗糙度缝隙内的气体密度也比主流区域的密度要大,这说明气体分子在粗糙壁面停留的时间较长,气体分子在粗糙壁面的集聚现象明显.图11为不同粗糙度高度下近壁面区域的正应力分布曲线,由于正应力的分布和气体密度直接相关,因此正应力分布曲线呈现出和密度分布相似的规律.

图10 不同壁面粗糙度下近壁区域流场密度分布Fig.10.Near-wall density distribution with different surface roughness.

图12为不同大小的粗糙度情形下剪切应力分布曲线.在主流区域剪切应力的大小趋于恒定,且随着粗糙度高度的增大而增大.当粗糙度高度增大到一定程度时,剪切应力不再继续增大,比如当粗糙度高度为0.54和0.81 nm时,剪切应力分别为18.85和18.99 kPa,二者几乎相等.通过剪切应力与切向动量适应系数的关系,计算得到气体分子对不同粗糙壁面的切向动量适应系数分别为0.87,0.98,0.99,可以看到当粗糙度大于0.54 nm时,TMAC的数值已经接近于1.0,表明气体分子的切向动量与粗糙壁面完全适应.上述结果说明粗糙度越大,气体与壁面的动量交换程度越高,气体越容易在壁面发生漫反射.

图11 不同壁面粗糙度下近壁区域流场正应力分布Fig.11.Near-wall normal stress distribution with different surface roughness.

图12 不同壁面粗糙度下近壁区域流场剪切应力分布Fig.12.Near-wall shear stress distribution with different surface roughness.

3.3 系统温度对流动的影响分析

为了研究系统温度对纳米通道内气体流动特性及切向动量适应系数的影响规律,针对通道高度为10.90 nm的气体剪切流动展开模拟,所有算例中的气体分子数目都保持一致,待系统平衡后,统计获得流场内速度、密度、正应力及剪切应力分布.

图13为不同温度条件下流场的无量纲速度分布,可知温度的变化明显地改变了流场的速度分布.为了更好地观察温度对流场速度的影响,截取了距离壁面2.0 nm的区域内的速度分布,如图14所示.从图14可以清晰地看到,整个流场的速度都随着温度的增大而减小.近壁面区域速度峰值大小的变化直观印证了这一点,由248 K时的51.54 m/s降低到了373 K时的41.03 m/s.此外,主流区域的速度剪切率也随着温度的增大而降低.

图13 不同温度条件下的流场无量纲速度分布Fig.13.Normalized velocity profiles with different temperature conditions.

图14 不同温度条件下近壁区域流场速度分布Fig.14.Near-wall velocity profiles with different temperature conditions.

图15为不同温度条件下的流场密度分布,由于主流区域的密度趋于定值,本文只展示了距壁面2.0 nm范围内的密度分布,后面的正应力分布及剪切应力也采用了这种表示方法.由图15可知,在壁面作用力影响的近壁面区域,气体密度随着温度的降低而升高,由于在初始设置计算条件时,所有算例都布置了相同数目的气体分子,因此,在温度较低时,气体分子在壁面的聚集会导致流场中心区域的密度有所降低.图16为流场中的正应力分布,正应力的分布主要和气体温度相关,流场的温度越高,则气体的正应力越大,MD的计算结果也呈现了这种趋势,即流场的正应力随着温度的升高而升高.将MD计算的主流区域的正应力结果和理想气体理论预测的结果列于表3,通过对比发现,对于所有温度条件下MD计算的结果和理论预测的误差不超过2%,以上结果也直接验证了MD在计算宏观量时的准确性.

图15 不同温度条件下近壁区域流场密度分布Fig.15.Near-wall density distribution with different temperature conditions.

图16 不同温度条件下近壁区域流场正应力分布Fig.16.Near-wall normal stress distribution with different temperature conditions.

表3 不同温度条件下主流区域正应力与理论值对比Table 3.Comparisons of normal stress with its theoretic values in bulk flow under different temperature conditions.

图17为不同温度条件下的剪切应力的分布,流场内的剪切应力在总体上是随着温度的升高而增大的.由气体动理论可知,剪切应力的变化不仅受到系统温度的影响,它还受到气体分子对壁面切向动量适应系数的影响[39].计算得到的切向动量适应系数如图18所示,随着温度的升高,气体分子对壁面的动量适应随之降低.温度除了对TMAC有影响,对气体分子在纳米通道内有效黏性系数也有明显的影响,由气体分子运动论可知,气体温度越高,则气体分子间的碰撞越频繁,则气体的黏性越大.气体的有效黏性系数可以通过下式计算：

图17 不同温度条件下近壁区域流场剪切应力分布Fig.17.Near-wall normal stress distribution with different temperature conditions.

图18 TMAC随温度的变化Fig.18.TMAC values with different temperature conditions.

其中τshear为流场中心区域剪切应力,∂u/∂z为流场中心区域的速度剪切率.根据(5)式计算得到了不同温度下气体有效黏性系数如图19所示,气体分子的有效黏性系数随着温度的升高而增大.

图19 有效黏性系数随温度的变化Fig.19.Effective viscosity values at different temperature.

4 结 论

采用三维非平衡分子动力学模拟方法研究了过渡区纳米通道内的气体剪切流动,获得了通道内气体的各物理量的分布,深入研究了近壁区域的流动特性,探讨了壁面条件及系统温度对纳米尺度气体流动的影响规律,并进一步分析了切向动量适应系数的变化趋势,发现了以下规律：

1)壁面原子的振幅相对于气体分子平均自由程来说是小量,壁面弹性系数的变化对流动的影响不明显,但是距壁面约一个分子直径处的速度峰值以及切向动量适应系数会随着壁面弹性系数的增大而减小;

2)壁面粗糙度对流动的影响主要通过切向动量适应系数的变化来体现,随着粗糙度的增大,气体分子对壁面的切向动量适应程度逐渐增强,直至完全适应;由于壁面粗糙度的存在会使壁面力场发生改变,进而造成近壁面区域速度、密度及应力分布的波动;

3)系统温度的变化对流动的影响是全局性的,无论是在近壁面区域还是在主流区域,气体的速度、密度、正应力和剪切应力都随着温度变化而改变,温度变化还会引起宏观参数的改变,气体分子的动量适应系数随着温度的升高而降低,有效黏性系数随着温度的升高而增大;

4)对于纳米尺度条件下的气体流动,由于表面积/体积之比极大,壁面力场对流场的主导作用更加明显,流动物理量对于壁面条件和流场温度的改变也十分敏感;壁面刚度的改变虽然对主流区域没有影响,但是对近壁区域的流动仍然有一定程度的干扰;而壁面粗糙度无论对主流区域还是近壁区域的流动都有较大影响,且气体分子在壁面处的聚集现象明显;系统温度的影响则是全局性的,可以显著改变气体与气体、气体与壁面之间的动量交换程度,因此对主流区域和近壁面区域的影响也都较大.