1992年, Qian等^[43]提出的LBM基本模型— — DnQm模型得到广泛应用。模拟非常规天然气流动时, 常采用BGK近似碰撞项^[44]的单松弛格子Boltzmann模型（LBGK）。LBGK模型是连续Boltzmann-BGK方程的特殊离散形式, 而连续Boltzmann-BGK方程是连续Boltzmann方程的简化。对于一个完整的LBGK模型, 应该包括演化方程、平衡态分布函数和格子模型, 其中演化方程表示如下：

${{f}_{i}}\left( X+{{c}_{i}}\delta t, t+\delta t \right)={{f}_{i}}\left( X, t \right)-$ $\frac{1}{\tau }\left[ {{f}_{i}}\left( X, t \right)-f_{i, eq}^{{}}\left( X, t \right) \right]+\delta t{{F}_{i}}$ （1）

在模拟二维空间流动时, 常采用九速离散的D2Q9格子模型, 此时, 平衡态分布函数为：

$f_{i, eq}^{{}}=\rho {{\omega }_{i}}\left[ 1+\frac{{{c}_{i}}\cdot u}{c_{s}^{2}}+\frac{{{\left( {{c}_{i}}\cdot u \right)}^{2}}}{2c_{s}^{4}}-\frac{{{u}^{2}}}{2c_{s}^{2}} \right]$ （2）

其中 ${{c}_{\text{s}}}=\frac{c}{\sqrt{3}}$ $c=\frac{\delta x}{\delta t}$ $\rho =\sum\limits_{i=0}^{8}{{{f}_{i}}}$

格子空间动量表示为：

$\rho u=\sum\limits_{i=0}^{8}{{{f}_{i}}{{c}_{i}}}+\frac{\delta t}{2}F$ （3）

在D2Q9格子模型中, 权重系数${{\omega }_{i}}$和离散格子速度c分别表示为：

${{\omega }_{i}}=\left\{ \begin{array}{* {35}{l}} \frac{4}{9} & c_{i}^{2}=0 \\ \frac{1}{9} & c_{i}^{2}={{c}^{2}} \\ \frac{1}{36} & c_{i}^{2}=2{{c}^{2}} \\ \end{array} \right.$ （4）

$c=c\left[ \begin{matrix} 0 & 1 & 0 & -1 & 0 & 1 & -1 & -1 & 1 \\ 0 & 0 & 1 & 0 & -1 & 1 & 1 & -1 & -1 \\ \end{matrix} \right]$ （5）

此外, 格子空间中运动黏度与无因次松弛时间满足关系式^[45]：

$\nu =c_{s}^{2}\left( \tau -\frac{1}{2} \right)\delta t$ （6）

上述黏度处理方式使得LBGK模型在计算时具有二阶精度。格子空间中, 单组分单相气体压力表示为^[28]：

$p=c_{s}^{2}\rho $ （7）

外加作用力的离散表达式^[28]：

${{F}_{i}}={{\omega }_{i}}\left( 1-\frac{1}{2\tau } \right)\left[ \frac{F\cdot {{c}_{i}}}{c_{s}^{2}}+\frac{\left( uF+Fu \right):\left( {{c}_{i}}{{c}_{i}}-c_{s}^{2}I \right)}{2c_{s}^{4}} \right]$ （8）

在宏观液体流动模拟中, 通常考虑雷诺数为流动特征控制参数, 可以通过雷诺数确定无因次松弛时间; 在微尺度气体流动模拟中, 努森数为流动特征控制参数, 此时需要对无因次松弛时间进行修正或重新定义。前人在这方面做了大量的研究, 通过建立努森数和无因次松弛时间的关系式（见表1）, 进而确定无因次松弛时间。2002年, Nie等^[30]和Lim等^[15]推导出无因次松弛时间表达式, 但二者存在明显差别, 前者有参数“ 1/2” 项, 说明模型具有二阶精度^[28], 而后者没有。Niu等^[14]推导出的无因次松弛时间表达式中, 当模拟气体为单原子气体时比热容为5/3, 此时τ ≈ Kn H, 显然该方法不具备二阶计算精度。根据郭照立等^[28]的推论, 无因次松弛时间可以表示为$\tau ={1}/{2}\; +\left( c/\bar{c} \right)Kn\text{ }N$, 不同文献中对于粒子平均速度$\bar{c}$有不同的选择, 从而导致无因次松弛时间存在不同表达形式。Lee等^[31]推导的公式中选取粒子平均速度近似于格子速度, 但是根据Guo等^[32]的研究, 该处理方式不满足相容性条件。Tang等^[33]和Zhang等^[34]所得公式中选取粒子平均速度为$\bar{c}=$$\sqrt{{8RT}/{\pi }\; }$。Guo等^{[32, 35]}从分子动力学理论出发, 近似推导出粒子平均速度表达式为$\sqrt{{\pi RT}/{2}\; }$, 从而得到无因次松弛时间表达式。目前, 研究者一般直接应用Tang等^[33]、Zhang等^[34]或Guo等^{[32, 35]}推导出的表达式或在此基础上进行改进。如在高温高压条件下, Shi等^[36]考虑气体黏度受稠密性影响对其做出相应修正, 此时可将Tang等^[33]、Zhang等^[34]和Guo等^[35]提出的无因次松弛时间表达式修正为$\tau ={1}/{2}\; +$$\sqrt{{3\pi }/{8}\; }Kn\text{ }N{{{\left[ 1+{2b{{\rho }_{r}}\chi }/{\left( D+2 \right)}\; \right]}^{2}}}/{\chi }\; $或$\tau ={1}/{2}\; +$ $\sqrt{{6}/{\pi }\; }Kn\text{ }N{{{\left[ 1+{2b{{\rho }_{r}}\chi }/{\left( D+2 \right)}\; \right]}^{2}}}/{\chi }\; $, 其中$\chi =1+$ $\left( {5}/{8}\; \right)b{{\rho }_{r}}+0.286\ 9{{\left( b{{\rho }_{r}} \right)}^{2}}+0.110\ 3{{\left( b{{\rho }_{r}} \right)}^{3}}+0.038\ 6{{\left( b{{\rho }_{r}} \right)}^{4}}$, $b={2\pi d_{r}^{3}}/{3{{m}_{r}}}\; $。

表1

表1

表1 文献中无因次松弛时间和努森数关系 文献 公式 Nie等[30] $\tau ={1}/{2}\; +Kn\text{ }H{\rho }/{a}\; $ Lim等[15] $\tau =Kn\text{ }N$ Niu等[14] $\tau =\sqrt{{6}/{\gamma \pi }\; }Kn\text{ }H$ Lee等[31] $\tau ={1}/{2}\; +Kn\text{ }N$ Tang等[33], Zhang等[34] $\tau ={1}/{2}\; +\sqrt{{3\pi }/{8}\; }Kn\text{ }N$ Guo等[35] $\tau ={1}/{2}\; +\sqrt{{6}/{\pi }\; }Kn\text{ }N$ Guo等[32] （努森层修正） $\tau ={1}/{2}\; +\sqrt{{6}/{\pi }\; }Kn\text{ }N\text{ }\psi \left( Kn \right)$, $\psi \left( Kn \right)=\left( {2}/{\pi }\; \right)\arctan \left( \sqrt{2}K{{n}^{{-3}/{4}\; }} \right)$ Li等[38] （努森层修正） $\tau ={1}/{2}\; +\sqrt{{6}/{\pi }\; }Kn\text{ }N\text{ }\psi \left( Kn \right)$, $\psi \left( Kn \right)={1}/{\left( 1+2Kn \right)}\; $ Zhang等[40] （努森层修正） $\tau ={1}/{2}\; +{\sqrt{{3\pi }/{8}\; }Kn\text{ }N}/{\left\{ 1+0.7\left[ {{\text{e}}^{{-Cy}/{\lambda }\; }}+{{\text{e}}^{{-C\left( H-y \right)}/{\lambda }\; }} \right] \right\}}\; $ Tang等[41] （努森层修正） $\tau ={1}/{2}\; +\left( {\lambda }/{{{\lambda }_{0}}}\; \right)\sqrt{{\pi }/{8}\; }\left( {c}/{{{c}_{s}}}\; \right)Kn\text{ }N$ Guo等[42] （努森层修正） $\tau ={1}/{2}\; +\sqrt{{6}/{\pi }\; }Kn\text{ }N\left[ 1+\left( {y}/{\lambda -1}\; \right){{\text{e}}^{{-y}/{\lambda }\; }}-{{\left( {y}/{\lambda }\; \right)}^{2}}Ei\left( {y}/{\lambda }\; \right) \right]$

表1 文献中无因次松弛时间和努森数关系

在微尺度流动模拟中, 还需要考虑边界努森层带来的影响。固体壁面的存在使得靠近壁面位置的分子平均自由程被截断, 导致一定范围内气体分子（即努森层）的运动规律受到影响^{[8, 11]}。当流动通道特征长度较大时, 努森层在整个流动通道中所占比例小, 其对流动的影响可以忽略不计; 但随着特征长度的减小, 努森层在流动通道内所占比例逐渐增大, 其对流动产生的影响也越来越大^[8], 产生微尺度流动特征。因此, 研究者对无因次松弛时间进行了相应的修正。总体来说, 主要通过修正分子平均自由程以表征边界努森层的影响。一些学者将努森层中气体分子平均自由程减小等价转换为整个流场中气体分子平均自由程的减小。Stops等^[37]提出了二维平板流动下的分子平均自由程修正式${{\lambda }_{e}}=\lambda \psi \left( Kn \right)$, 但修正系数$\psi \left( Kn \right)$过于复杂, 实际计算不方便; Guo等^[32]在此基础上提出新的修正系数表达式; Li等^[38]根据气体动力黏度与分子平均自由程的正比例关系, 采用Michalis等^[39]提出的Bosanquet-type型修正系数。Zhang等^[40]、Tang等^[41]、Guo等^[42]考虑流场中不同位置与壁面的距离对该位置处分子平均自由程进行修正, 推导出相应的无因次松弛时间。

上述无因次松弛时间表达式均通过分子动力学理论推导而来, 在微尺度流动模拟研究中应用广泛。如图1所示, 作出部分表达式中无因次松弛时间随努森数的变化曲线, 可以看到在相同网格数及努森数条件下, 无因次松弛时间值存在差异, 并且在过渡流区尤其明显。无因次松弛时间对微尺度气体流动模拟研究影响重大, 微小的取值差别会导致模拟结果出现差异, 可能造成模拟结果偏离正确值。

图1

	Figure Option ViewDownloadNew Window
	图1 不同努森数下无因次松弛时间取值（N=100）

前人在推导无因次松弛时间时, 分子平均自由程与松弛时间、粒子平均运动速度满足关系式^[36]：

$\lambda ={{\tau }_{\text{s}}}\bar{c}$ （9）

动力黏度和压力、松弛时间满足关系^[32]：

$\mu =p{{\tau }_{\text{s}}}$ （10）

考虑高温高压非理想气体稠密性修正时^[36], 动力黏度可表示为：

$\mu =p\frac{\lambda }{{\bar{c}}}{{{\left[ 1+{2b\rho \chi }/{\left( D+2 \right)}\; \right]}^{2}}}/{\chi }\; $ （11）

不考虑壁面产生的影响, 绘制出温度为350 K时（11）式对应的真实物理空间中甲烷气体分子动力黏度与压力（1~100 MPa）关系图（见图2）, 并与美国国家标准与技术研究院（NIST）化学数据库中的数据对比。从图中可以看出, 压力高于20 MPa时, 各表达式计算的甲烷气体分子动力黏度均与NIST数据库中数据存在较大偏差。此时, 上述无因次松弛时间表达式无法反映气体真实情况。

图2

	Figure Option ViewDownloadNew Window
	图2 不同表达式计算的的甲烷动力黏度随压力的变化

在微尺度流动模拟中, 应使模拟的格子空间量和真实物理空间量满足相似准则^{[46, 47]}, 如马赫数、努森数及雷诺数相等, 才能真正反映实际情况, 揭示客观规律。在宏观液体流动模拟中, 考虑数值模拟条件与真实流动条件下的雷诺数相等, 实际上是满足格子空间和真实物理空间中的惯性力与摩擦力之比相等。值得注意的是, 对于宏观大尺度液体流动, 可忽略液体的压缩性, 但针对微观气体流动应该考虑气体的压缩性。此时, 需要严格满足模拟条件下的马赫数和实际马赫数相等, 也就是满足二者的惯性力和弹性力之比相等。

由格子空间和真实物理空间下的雷诺数相等可得：

$\frac{uH}{\nu }=\frac{{{u}_{r}}{{H}_{r}}}{{{\nu }_{r}}}$ （12）

格子空间和真实物理空间中的努森数相等, 则：

$\frac{\lambda }{H}=\frac{{{\lambda }_{r}}}{{{H}_{r}}}$ （13）

满足格子空间和真实物理空间中的马赫数相等, 则：

$\frac{u}{{{c}_{s}}}=\frac{{{u}_{r}}}{{{c}_{sr}}}$ （14）

联立（6）式及（12）— （14）式可得：

$\tau =\frac{1}{2}+\frac{\lambda {{\nu }_{r}}}{{{c}_{s}}{{c}_{sr}}{{\lambda }_{\text{r}}}\delta t}$ （15）

将$\lambda =Kn\text{ }H$、$H=N\delta x$和${{c}_{s}}=\frac{c}{\sqrt{3}}$带入（15）式可得：

$\tau =\frac{1}{2}+\frac{\sqrt{3}{{\nu }_{r}}}{{{\lambda }_{r}}{{c}_{sr}}}Kn\text{ }N$ （16）

真实物理空间中气体分子平均自由程计算式为：

${{\lambda }_{r}}=\frac{{{m}_{r}}}{\sqrt{2}\pi {{\rho }_{r}}d_{r}^{2}}$ （17）

前人的研究结果表明, 可以通过修正无因次松弛时间表达式以表征微尺度条件下固体壁面附近努森层对流动的影响^{[32, 38]}：

$\tau =\frac{1}{2}+\frac{\sqrt{3}{{\nu }_{r}}}{{{\lambda }_{r}}{{c}_{sr}}}Kn\text{ }N\text{ }\psi \left( Kn \right)$ （18）

为便于模拟计算, 对于上式中修正系数$\psi \left( Kn \right)$, 采用Guo等^[32]提出的$\psi \left( Kn \right)=\left( {2}/{\pi }\; \right)\arctan \left( \sqrt{2}K{{n}^{{-3}/{4}\; }} \right)$或Michalis等^[39]提出的$\psi \left( Kn \right)={1}/{\left( 1+2Kn \right)}\; $, 分别记为ψ （Kn）_1和ψ （Kn）_2。计算出上述两项修正系数在Kn为1× 10^-5~10时的取值分布, 可以看出二者取值存在明显差异（见图3）。在连续流区域, 两项修正系数基本保持不变, 并且存在ψ （Kn）_1≈ ψ （Kn）_ 2≈ 1.0; 在滑移流和过渡流区域, 修正系数开始下降, 且下降速率先增加后减小。修正系数间差异主要出现在滑移流区和过渡流区, 同时在这两个区间内ψ （Kn）_1> ψ （Kn）_2。此外, 在Kn≈ 1时存在最大差值Δ ψ （Kn）_max, 其中$\Delta \psi \left( Kn \right)\text{=}\psi \left( Kn \right)\_1-\psi \left( Kn \right)\_2$。由于（18）式中的气体运动黏度可以根据温度和压力从NIST数据库内查询, 无需再对气体稠密性影响进行单独修正。由于无因次松弛时间的微小差别会导致模拟结果产生明显差异, 关系到流动模拟能否正确捕捉微尺度下的物理现象, 因此需要优选修正系数。

图3

	Figure Option ViewDownloadNew Window
	图3 不同努森数下修正系数的取值及差异

一般微尺度模拟中, 考虑气体在微通道内周期流动时, 通常在流动方向进出口处采用周期边界条件^[47]; 考虑气体在微通道中受进出口压差驱动流动时, 在进出口处采用Zou等^[48]提出的压力边界条件或Guo等^[49]提出的非平衡态外推边界条件。相比之下孔道上下壁面处边界处理更复杂, 通过在壁面处采用合理的边界条件可以有效捕捉边界滑移等现象。

目前, 处理壁面处微尺度气体流动时通常采用两种边界条件^{[7, 35]}：组合反弹/镜面反射边界条件（CBBSR）和离散Maxwell边界条件（DM）。这两种边界条件在本质上是等价的, 本文采用CBBSR边界条件^[35]。以二维微通道下边界为例, 网格划分如图4所示, 本文将网格节点设置在壁面边界位置处有利于直接捕捉边界滑移速度。经过碰撞步和迁移步后, 未知分布函数为${{f}_{2}}$、${{f}_{5}}$和${{f}_{6}}$, 采用CBBSR边界条件可以得到：

$\left\{ \begin{matrix} {{f}_{2}}={{{\bar{f}}}_{4}} \\ {{f}_{5}}=r{{{\bar{f}}}_{7}}+\left( 1-r \right){{{\bar{f}}}_{8}} \\ {{f}_{6}}=r{{{\bar{f}}}_{8}}+\left( 1-r \right){{{\bar{f}}}_{7}} \\ \end{matrix} \right.$ （19）

图4

	Figure Option ViewDownloadNew Window
	图4 网格划分示意图（数字表示离散速度方向）

在微尺度流动模拟中, 组合系数r的取值直接影响模拟得到的边界滑移速度, 若将r设为定值^[50]则无法反映出模拟条件的差异。为消除数值离散效应并使模拟满足二阶滑移边界条件, 根据Guo等^[35]的相关理论, 结合本文提出的无因次松弛时间推导得出组合系数表达式：

$r=\left\{ 1-\frac{4{{\nu }_{r}}}{\sqrt{3}{{\lambda }_{r}}{{c}_{sr}}}Kn\text{ }\psi \left( Kn \right)+ \right.$ ${{\left. \frac{{{\lambda }_{r}}{{c}_{sr}}}{\sqrt{3}{{\nu }_{r}}}\left[ {{A}_{1}}+2{{A}_{2}}Kn\text{ }\psi \left( Kn \right)+\frac{1}{4{{N}^{2}}Kn\text{ }\psi \left( Kn \right)} \right] \right\}}^{-1}}$ （20）

Guo等^[51]在广义二阶滑移边界条件中给出：

${{A}_{1}}=\frac{2-\sigma }{\sigma }\left( 1-0.181\text{ }7\text{ }\sigma \right)$ ${{A}_{2}}={{\sigma }^{2}}\left( \frac{1}{\text{ }\!\!\pi\!\!\text{ }}+\frac{1}{2}A_{1}^{2} \right)$

其中 $\sigma ={\left( 2-{{\sigma }_{\text{v}}} \right)}/{{{\sigma }_{\text{v}}}}\; $

首先, 利用等温条件下甲烷气体在无限长微通道中的体积力驱动流动进行模型验证。通过在进出口设置周期流动边界条件^[47], 可以将气体流动等效为无限长微通道中的流动。通常将压力梯度设置为很小的常量^[38], 并且将整个流场中外力项统一为$F=\left( {{F}_{x}}, {{F}_{y}} \right)$, 本模拟取F_x为10^-6, F_y为零。此时, 沿流动方向气体压力变化微小, 可以忽略气体宏观物性参数差异。将整个流场中的无因次松弛时间设为恒定, 同时可以忽略压缩效应的影响, 只需要考虑气体的滑脱效应和边界努森层的影响。网格划分为N_x× N_y=30× 300, 由于采用周期流动边界条件, 流动方向（即x方向）网格数少一些。（18）式和（20）式中涉及到的所有甲烷气体物性参数均来自NIST标准参考数据库, 通过改变模拟的微通道上下壁面间距H_r可以得到不同Kn及对应的无因次松弛时间和组合系数。由于需要模拟较高Kn下的流动, 为保证通道特征长度为纳米级, 模拟压力不能太高, 因此取T_r为400 K, p_r为0.1 MPa。

Ohwada等^[52]通过直接求解线性化Boltzmann方程对二维微通道中气体进行流动模拟, 通常将其结果用于验证微尺度流动模型的准确性^{[38, 40, 51, 53]}。采用本文模型进行模拟计算, 得到不同Kn条件下出口处的无因次速度剖面U/U_ave, 其中${{U}_{\text{ave}}}=\left( {1}/{H}\; \right)\int_{0}^{H}{Udy}$。采用不同的修正系数ψ （Kn）_1和ψ （Kn）_2时, 相应的模拟计算结果分别记为“ Present LB_1” 和“ Present LB_2” 。将模拟结果同Ohwada等的计算结果进行对比以验证模型准确性并优选修正系数, 同时还将计算结果同Zhang等^[40]、Li等^[38]和Guo等^[53]的计算结果进行对比（见图5）。

图5

	Figure Option ViewDownloadNew Window
	图5 甲烷气体在微通道中周期流动时出口处无因次速度剖面

由图5可以看出, 当$\psi \left( Kn \right)=\psi \left( Kn \right)\_1$时, Kn取值不同本文模型计算精度不同。当Kn值为0.112 8时, 本文模型计算结果与Ohwada等直接求解线性化Boltzmann方程得到的结果具有较好的一致性, 但是在中心线位置处本文模型计算结果大于Zhang等和Guo等的结果; 随着Kn值的增大, 本文模型计算结果逐渐偏离Ohwada等的结果, 具体表现为边界处滑移速度偏大, 而通道中间速度偏小, 与Zhang等的计算结果呈现出很好的一致性; 当$Kn> 2.5$, 模拟不收敛。当$\psi \left( Kn \right)=\psi \left( Kn \right)\_2$时, 本文模型计算结果与Ohwada等的模拟结果始终保持很好的一致性, 说明此时模型具有更高的计算精度, 可以在较大的Kn变化范围内有效捕捉微尺度气体流动边界滑移现象。此外, 本文模拟结果与Li等（Li等使用了相同的修正系数$\psi \left( Kn \right)\_2$）的计算结果吻合度较高, 说明$\psi \left( Kn \right)=\psi \left( Kn \right)\_1$对边界滑移速度的准确捕捉具有重要影响。本文模型计算结果与Guo等计算结果始终存在偏差, 并且在Kn值较大时偏差非常明显。可以发现, 随着Kn值的增加, 边界滑移速度显著增大, 滑脱效应对流动的影响越来越明显。

为进一步验证本文模型的正确性并优选修正系数, 计算出不同Kn值条件下微通道内无因次流量^[38]$Q={\int_{0}^{H}{udy}}/{\left( {{{F}_{x}}{{H}^{2}}\sqrt{{RT}/{2}\; }}/{p}\; \right)}\; $, 并将计算结果同文献中数据进行对比（见图6）。可以看出, 当$\psi \left( Kn \right)=$ $\psi \left( Kn \right)\_1$, $Kn< 1.0$时, 本文模型计算结果与Ohwada等^[52]、Cercignani等^[54]、Li等^[38]的模拟结果以及真实气体的实验数据^[53]具有很好的一致性; 但随着Kn值增大, 本文模型计算结果逐渐增大并偏离上述文献中的模拟结果, 且当$Kn> 2.5$后计算不收敛, 无法模拟。当$\psi \left( Kn \right)=\psi \left( Kn \right)\_2$时, 在整个模拟的Kn取值范围内, 本文计算结果与上述文献中的模拟结果及实验数据具有较好的一致性, 从而验证了模型的正确性。通过对不同修正系数下的模拟计算结果进行对比, 表明采用修正系数$\psi \left( Kn \right)=\psi \left( Kn \right)\_2$可以使本文模型具有更高的计算精度和更广的适用范围。

图6

	Figure Option ViewDownloadNew Window
	图6 甲烷在微通道中周期流动时无因次流量随努森数的变化

在上述两种修正系数下计算出不同Kn值条件下出口处的归一化边界滑移速度${{{U}_{\text{slip}}}}/{{{U}_{0}}}\; $。将计算结果同Arkilic等^[55]考虑一阶滑移边界推导得到的解析解以及Nie等^[30]的结果进行对比（见图7）。可以看出, 在滑移流区内采用不同修正系数得到的计算结果基本吻合, 但随着Kn值的增加二者差距逐渐显现; 在过渡流区内, 二者差异进一步增加。同时, 在整个Kn取值范围内, 当$\psi \left( Kn \right)=\psi \left( Kn \right)\_1$时, 本文模型计算结果同Arkilic等的解析解具有更好的一致性。实际上, 在Arkilic等的推导过程中仅考虑一阶滑移边界条件, 因此采用$\psi \left( Kn \right)\_1$时本文模型只具有一阶计算精度。在Nie等的研究中, 采用的是简单的反弹边界条件, 可以看到在滑移流区, Nie等的计算结果小于其余3种计算结果, 在过渡流区则相反。当$\psi \left( Kn \right)=\psi \left( Kn \right)\_2$时, 过渡流区本文模型模拟计算结果小于其余3种结果, 说明采用$\psi \left( Kn \right)\_1$导致计算得到的边界滑移速度在过渡流区偏大。

图7

	Figure Option ViewDownloadNew Window
	图7 甲烷在微尺度二维平板中周期流动时归一化边界滑移速度随努森数的变化

采用本文模型模拟甲烷在长直微通道中受进出口压差驱动的流动。对于微尺度气体流动, 在研究范围内温度压力变化微小, 气体流动速度缓慢, 属于典型的低马赫数流动。虽然理论上可以将微尺度气体流动视为不可压缩流动, 但是在整个宏观尺度上来看, 压缩性是客观存在的, 忽略压缩效应的模拟必然存在误差。在研究此类问题时, 需要同时考虑压缩效应、气体滑脱效应和边界努森层的影响。由于气体的可压缩性, 微通道内沿流动方向压力表现为非线性分布。将本文模拟结果同直接模拟蒙特卡洛方法（DSMC）计算结果、信息保存-直接模拟蒙特卡洛方法（IP-DSMC）计算结果、Li等^[38]、Guo等^[53]、Shen等^[56]的LBM模拟结果以及Arkilic等^[55]的解析解进行对比, 以验证模型的准确性。在上述文献中, 进出口压力比为1.4:1.0或2.0:1.0, 考虑模拟的压力梯度不能太大, 设定出口压力为0.01 MPa, 模拟温度恒定为400 K, 通过改变微通道上下壁面间距以实现不同Kn值条件下的流动模拟。气体沿流动通道方向压力发生变化, 此时在整个流场内气体的宏观物性参数也会发生相应变化。根据（18）式和（20）式, 可以看出不同位置处的无因次松弛时间、努森数以及组合系数不同, 在模拟时需要考虑各参数值随位置的变化。流场任意位置处努森数表示为$Kn\left( x, y \right)={K{{n}_{out}}{{\rho }_{out}}}/{\rho \left( x, y \right)}\; $。值得注意的是, 在计算无因次松弛时间和组合系数时涉及到${{{\nu }_{\text{r}}}}/{{{c}_{\text{sr}}}}\; $, 需要进行相应的预处理。通过NIST数据库查询得到${{{\nu }_{\text{r}}}}/{{{c}_{\text{sr}}}}\; $和${{\rho }_{\text{r}}}$在进出口压力范围内的一系列离散值, 然后拟合出二者关系式${{{\nu }_{\text{r}}}}/{{{c}_{\text{sr}}}}\; $=3× 10^-8/ρ _r。实际模拟时, 设置格子空间密度和真实物理空间密度在数值上相等, 此时可结合上述关系式计算出任意位置处的无因次松弛时间和组合系数。

取$\psi \left( Kn \right)=\psi \left( Kn \right)\_2$, 采用本文模型计算出$K{{n}_{\text{out}}}$分别为0.019 4、0.194 0和0.388 0时的沿程无因次压力偏移量${\left( p-{{p}_{line}} \right)}/{{{p}_{out}}}\; $。网格划分为N_x N_y=2 000× 20, 与参考文献[55]一致。

从图8可以看出, 采用本文模型计算得到的结果明显优于Arkilic等、Shen等和Guo等的结果, 与Li等的计算结果相近; 但相比之下本文模型计算结果与DSMC、IP-DSMC法的模拟结果更具有一致性, 说明本文模型能更为准确地捕捉微尺度气体的压缩效应。

图8

	Figure Option ViewDownloadNew Window
	图8 本文模型计算的长直微通道中气体沿程无因次压力偏移量与文献数据对比