1、1微尺度下气体在过渡区内流动的格子 Boltzmann 模拟摘要： 为研究微尺度下气体在过渡区内的流动特性，基于气体动理学及 Knudsen层效应理论，推导了 Knudsen数与无量纲松弛时间的关系；应用 Succi的边界处理方法和广义二阶滑移边界条件，推导了壁面滑移速度和反弹比例系数的计算公式，建立了适用于过渡区微尺度气体流动的格子 Boltzmann模型，并应用该模型对过渡区内微尺度Poiseuille流动进行模拟.结果表明，当稀薄参数取 1.64时，计算得到的无量纲速度剖面在整个过渡区与 Karniadakis给出的无量纲速度剖面吻合较好，无量纲速度分布在过渡区基本上保持为抛物线形状，边
2、界上的无量纲滑移速度随着 Knudsen数的增加而增大，中心线上的无量纲速度随着 Knudsen数的增加而减小. 关键词： 微尺度气体流动；格子 Boltzmann模型；Knudsen 数；滑移速度；过渡区；稀薄参数 中图分类号： V211.25文献标志码： ALattice Boltzmann Simulation of Microscale Gas Flows in the Transitional Regime LIU Jiali，ZHANG Jiye，ZHANG Weihua （Traction Power State Key Laboratory， Southwest Jiaoton
3、g University， Chengdu 610031， China） Abstract：In order to study the flow characteristics of 2microscale gas in the transitional regime， the relationship between Knudsen number and dimensionless relaxation time was derived based on the gas kinetic theory and the effect of Knudsen layer. Computational
4、 formulas for the slip velocity on the wall and the bounceback fraction were derived under a generalized secondorder slip boundary condition using the boundary treatment method proposed by Succi. Then， a lattice Boltzmann model for microscale gas flows in the transitional regime was established， and
5、 the microscale Poiseuille flows in the transitional regime were simulated. Computational results show that when the rarefaction parameter is equal to 1.64， the computed dimensionless velocity profile is in good agreement with the dimensionless velocity profile given by Karniadakis in the whole tran
6、sitional regime. The dimensionless velocity profile remains essentially a parabolic shape in the transitional regime. As Knudsen number increases， the dimensionless slip velocity rises in the boundary and falls in the center line. Key words：microscale gas flow； lattice Boltzmann model； Knudsen numbe
7、r； slip velocity； transitional regime； rarefaction parameter 微型化是自然科学和工程技术发展的重要趋势，尤其是微机电系统3技术的飞速发展，推动了这一研究成为热潮.发展关于微器件的基础科学和工程技术是蕴含在这些新技术中的需求，微机电系统中气体的特性是其中非常关键的问题，也是其它相关技术的基础.微机电系统中气体的流动会出现与常规大尺度气体流动1明显不同的现象，如边界滑移.微尺度气体流动的奇异特性主要由气体稀薄性引起，根据气体分子动力学，稀薄气体分子的速度分布函数满足 Boltzmann方程，其求解方法主要有分析方法和数值方法两大类.分析方
8、法主要包括线化 Boltzmann方程方法2、矩方法3和模型方程方法4.数值方法主要有速度滑移的NavierStokes方程求解56、直接模拟 Monte Carlo（direct simulation Monte Carlo， DSMC）方法7和基于 DSMC的信息保存法89. 近年来，格子 Boltzmann方法也被用于微尺度气体流动的模拟，这方面的研究始于 2002年 Nie10和 Lim11的研究工作.然而， Shen12指出当 Knudsen数较大（Kn0.1）时， Nie的模型会存在一个较大的压力负偏差. 随着模型的改进和发展，格子 Boltzmann方法在微尺度气体流动上的应用也
9、逐步完善1314.在 2006年以后，格子 Boltzmann方法开始用于过渡区流动的模拟1517，但也仅是在 Kn 本文在标准格子Boltzmann模型的基础上，从气体动理学理论和 Knudsen效应出发，推导出适用于过渡区微尺度气体流动的无量纲松弛时间与 Knudsen数的关系，进而将二阶速度滑移边界条件推广为广义二阶速度滑移边界条件，并基于 Succi的边界处理方法，推导出壁面滑移速度和反弹比例系数的计算4公式，建立了适用于过渡区微尺度气体流动的格子 Boltzmann模型；以过渡区的微尺度 Poiseuille流动为例进行数值分析，以验证本文建立的计算模型的正确性.1 过渡区微尺度气体
10、流动的格子 Boltzmann模型 1.1格子 Boltzmann模型格子 Boltzmann方程可以从描述稀薄气体流动的连续 Boltzmann方程得到，其一般形式为 f（x+et，t+t）-f（x，t）= -1f（x，t）-f，eq（x，t）+tF（x，t） ， （1） 式中： f（x，t）表示 t时刻在 x处的流体粒子在 方向上的速度分布函数，可简记为 f； f，eq（x，t）表示 t时刻在 x处的流体粒子在 方向上的平衡态分布函数，可简记为 f，eq； F（x，t）表示 t时刻在 x处的流体粒子在 方向上的外力，可简记为 F； 表示无量纲的松弛时间； t 表示时间步长. 离散速度模型采
11、用文献19中提出的 DdQm系列模型，其中， d表示空间维数， m 表示离散速度个数. 对于二维问题，一般采用 D2Q9模型，其离散速度为 e=（0，0） ，=0， c（cos（-1）/2，sin（-1）/2） ，=1，2，3，4， 2（cos（2-1）/2，sin（2-1）/2） ，=5，6，7，8， （2） 5式中： e 表示流体粒子在 方向上的离散速度； c=x/t， 其中： x 和 t 分别表示网格步长和时间步长. 相应的平衡态分布函数为 f，eq=1+e，uc2s+（eu）22c4s-u22c4s， （3） 式中： cs=1/3 表示格子声速； 表示权系数， =0，1，2，8； 0=
12、4/9； i=1/9（i=1，2，3，4） ， i=1/36（i=5，6，7，8） ； 表示宏观的流体密度； u 表示宏观的流体速度. 外力项 F（x，t）的模型20为 F=1-12euc2s+（eu）ec4sa， （4） 式中： a 表示外力加速度. 宏观的流体密度和流体速度可通过下式计算， =f， u=1fe+12ta.（5）1.2 无量纲松弛时间与 Knudsen数的6关系从动理学理论可知，气体分子平均自由程为 =2/（） ， 其中： 表示动力粘性系数； =8RT/ 表示气体分子平均速度； R 表示气体常数； T 表示温度. 另一方面，气体分子的平均自由程还可以表示为 =. 在格子 Bo
13、ltzmann模型中，格子速度表示为 c=RT， 其中： 是与模型相关的常数，对于 D2Q9模型， =3. 此外，需要注意的是 Boltzmann方程离散时，需要对 进行修正，修正量为-0.5，由此可得无量纲松弛时间与 Knudsen数（Kn）的关系为 =12+64 Kn， （6） 式中： 表示计算网络， =1/L，其中， L 为流动特征长度. 对于微尺度气体流动，气体在流经固体壁面时，在壁面附近会产生一个厚度为分子平均自由程量级的 Knudsen边界层，在该区域内分子间的碰撞频率大大降低，分子运动远远没有达到局部热力学平衡状态，应力与应变率不再满足通常的线性关系.为使格子 Boltzmann
14、模型能够描述Knudsen层内的非线性速度分布，通过引入 Knudsen层的有效粘性系数7e 对模型进行修正，在 Knudsen层内， e 可表示为 e=， （7） 式中： 表示修正函数. Beskok 和 Karniadakis建议在过渡区采用依赖于 Knudsen数的修正函数，其表达式为21 （Kn）=11+rKn， （8） 式中： r 表示稀薄参数. Vasilis22采用 DSMC方法，对过渡区的 Couette流和 Poiseuille流进行数值模拟，指出当 r=1.5时，此修正函数的计算结果在过渡区与DSMC的计算结果吻合很好. Knudsen 层对稀薄气体流动的影响可通过修正函数
15、对粘性系数的修正，将式（6）引入到格子 Boltzmann模型中，以有效粘性系数 e 代替粘性系数 ，则可得考虑 Knudsen层效应后，无量纲松弛时间与 Knudsen数的关系为 =12+64 Kn（Kn）.（9）1.3 微尺度气体流动的边界条件对于微尺度气体流动，气体在边界上存在速度滑移现象，必须采用能反映速度滑移效应的边界格式.本文采用 Succi提出的混合边界条件23，即将无滑移的反弹格式与无穷滑移的镜面反弹格式组合起来的混合格式，简记为 BSR格式. 考虑如图 1所示的气体在二维通道中的流动，假设密度 为常数， 8y方向的速度分量为 0，并且/x=0（表示任一物理量） ，外力加速度
16、a与壁面平行，采用半步长反弹的 BSR格式，即壁面置于 j=1/2及 j=n+1/2处. 图 1二维通道中 D2Q9模型的边界布置 Fig.1Boundary setting of D2Q9 model in two dimensional channel 在时刻 t，首先执行碰撞过程 f （x，t）=f（x，t）+ 1f（x，t）-f，eq（x，t）+tF， （10） 式中： f（x，t）表示碰撞后的速度分布函数. 对于 1jN的格点可以执行迁移过程为 f（x+et，t+t）=f （x，t）.（11） 对于 j=1的格点，只有 f0、f1、f3、f4、f7 和 f8可以通过迁移过程得到， f
17、2、f5 和 f6需要根据边界条件确定； 对于 j=n的格点， f0、f1、f2、f3、f5 和 f6可以通过迁移过程得到， f4、f7 和 f8需要根据边界条件确定，采用 BSR格式，其表达式分别为 f2=f4， f5=rbf7+（1-rb）f8， f6=rbf8+（1-rb）f7， f4=f2， f7=rbf5+（1-rb）f6， 9f8=rbf6+（1-rb）f5， （12） 式中： rb0，1表示反弹比例系数，即粒子在与壁面作用时沿原路弹回所占的比例. 对于定常外力驱动下的 Poiseuille流动， BSR格式在边界上产生一个无量纲滑移速度（采用中心线速度进行无量纲化）24 Us=2
18、（1-rb）rb（2-1）+ 43（2-1）22-2， （13） 式中： Us 表示由边界条件确定的滑移速度， Us=us/uc， 其中， uc 表示中心线上的最大速度. 根据式（9） 、 （13） ，Us 还可采用 Knudsen数表示为 Us=1-rbrb6Kn（Kn）+ 2Kn22（Kn）-2.（14） 考虑的二阶滑移边界条件为 us=C1unwall-C22unwall， （15） 式中： C1、C2 表示依赖气体与壁面相互作用的参数. Knudsen 层对二阶滑移边界条件的影响可通过修正函数对粘性系数的修正，由式（7）引入到微尺度气体流动的边界条件模型中，广义二阶滑移边界条件为 10
19、us=C1（Kn）unwall-C222（Kn）unwall.（16） 利用广义二阶滑移边界条件，可得 Poiseuille流动在壁面上的无量纲滑移速度（采用中心线速度进行无量纲化）为 Us=4C1Kn（Kn）+8C2Kn22（Kn）.（17） 为实现上述滑移速度，由式（14） 、 （17）可知， BSR格式中的参数rb必须选择为 rb=1+162Kn（Kn）+4C1+ （8C2-2）Kn（Kn）-1， （18） 式中： C1=0.818 3； C2=0.653 1.2数值计算管道的稀薄气体流动广泛应用于各种技术，如真空管道、微机电系统、空天飞行器等.在管道两端以外力驱动的Poiseuille流动，引起了许多科学家和技术人员的研究兴趣.由于对这类流动的研究比较透彻，所以理论及试验研究结果可以作为检验新模型和新算法的参考基准.两个平行平板之间的微尺度 Poiseuille流动如图 1所示. Karniadakis采用局部平均速度对速度分布进行无量纲化，得到了仅依赖于 Knudsen数和 y的无量纲速度剖面的表达式25为 U（y，Kn）=-y2L2+yL+Kn1-bKn16+Kn1-bKn， （19） 式中： y 表示距通道下壁面的距离. 当 b=-1时，用式（19）可导出较为准确的适用于过渡区流动的速度分布.
