摘要： 為研究微尺度下氣體在過渡區(qū)內(nèi)的流動(dòng)特性，基于氣體動(dòng)理學(xué)及Knudsen層效應(yīng)理論，推導(dǎo)了Knudsen數(shù)與無量綱松弛時(shí)間的關(guān)系；應(yīng)用Succi的邊界處理方法和廣義二階滑移邊界條件，推導(dǎo)了壁面滑移速度和反彈比例系數(shù)的計(jì)算公式，建立了適用于過渡區(qū)微尺度氣體流動(dòng)的格子Boltzmann模型，并應(yīng)用該模型對(duì)過渡區(qū)內(nèi)微尺度Poiseuille流動(dòng)進(jìn)行模擬.結(jié)果表明，當(dāng)稀薄參數(shù)取1.64時(shí)，計(jì)算得到的無量綱速度剖面在整個(gè)過渡區(qū)與Karniadakis給出的無量綱速度剖面吻合較好，無量綱速度分布在過渡區(qū)基本上保持為拋物線形狀，邊界上的無量綱滑移速度隨著Knudsen數(shù)的增加而增大，中心線上的無量綱速度隨著Knudsen數(shù)的增加而減小.

關(guān)鍵詞： 微尺度氣體流動(dòng)；格子Boltzmann模型；Knudsen數(shù)；滑移速度；過渡區(qū)；稀薄參數(shù)

中圖分類號(hào)： V211.25文獻(xiàn)標(biāo)志碼： ALattice Boltzmann Simulation of Microscale

Gas Flows in the Transitional Regime LIU Jiali，ZHANG Jiye，ZHANG Weihua

（Traction Power State Key Laboratory， Southwest Jiaotong University， Chengdu 610031， China）

Abstract：In order to study the flow characteristics of microscale 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 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 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 transitional 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 number； slip velocity； transitional regime； rarefaction parameter

微型化是自然科學(xué)和工程技術(shù)發(fā)展的重要趨勢(shì)，尤其是微機(jī)電系統(tǒng)技術(shù)的飛速發(fā)展，推動(dòng)了這一研究成為熱潮.發(fā)展關(guān)于微器件的基礎(chǔ)科學(xué)和工程技術(shù)是蘊(yùn)含在這些新技術(shù)中的需求，微機(jī)電系統(tǒng)中氣體的特性是其中非常關(guān)鍵的問題，也是其它相關(guān)技術(shù)的基礎(chǔ).微機(jī)電系統(tǒng)中氣體的流動(dòng)會(huì)出現(xiàn)與常規(guī)大尺度氣體流動(dòng)[1]明顯不同的現(xiàn)象，如邊界滑移.微尺度氣體流動(dòng)的奇異特性主要由氣體稀薄性引起，根據(jù)氣體分子動(dòng)力學(xué)，稀薄氣體分子的速度分布函數(shù)滿足Boltzmann方程，其求解方法主要有分析方法和數(shù)值方法兩大類.分析方法主要包括線化Boltzmann方程方法[2]、矩方法[3]和模型方程方法[4].數(shù)值方法主要有速度滑移的NavierStokes方程求解[56]、直接模擬Monte Carlo（direct simulation Monte Carlo， DSMC）方法[7]和基于DSMC的信息保存法[89].

近年來，格子Boltzmann方法也被用于微尺度氣體流動(dòng)的模擬，這方面的研究始于2002年Nie[10]和Lim[11]的研究工作.然而， Shen[12]指出當(dāng)Knudsen數(shù)較大（Kn>0.1）時(shí)， Nie的模型會(huì)存在一個(gè)較大的壓力負(fù)偏差.

隨著模型的改進(jìn)和發(fā)展，格子Boltzmann方法在微尺度氣體流動(dòng)上的應(yīng)用也逐步完善[1314].在2006年以后，格子Boltzmann方法開始用于過渡區(qū)流動(dòng)的模擬[1517]，但也僅是在Kn<0.5的情況下與DSMC的結(jié)果吻合的比較好.文獻(xiàn)[18]中提出了包含多個(gè)有效松弛時(shí)間的MRTLBE模型，在更大的Knudsen數(shù)范圍內(nèi)仍能取得了較好的模擬結(jié)果，但模型相西南交通大學(xué)學(xué)報(bào)第48卷第4期劉加利等：微尺度下氣體在過渡區(qū)內(nèi)流動(dòng)的格子Boltzmann模擬對(duì)復(fù)雜.從理論上講，格子Boltzmann方法是由Boltzmann方程的簡化模型——BoltzmannBGK方程發(fā)展而來的，該方法適用于過渡區(qū)流動(dòng)的模擬.

本文在標(biāo)準(zhǔn)格子Boltzmann模型的基礎(chǔ)上，從氣體動(dòng)理學(xué)理論和Knudsen效應(yīng)出發(fā)，推導(dǎo)出適用于過渡區(qū)微尺度氣體流動(dòng)的無量綱松弛時(shí)間與Knudsen數(shù)的關(guān)系，進(jìn)而將二階速度滑移邊界條件推廣為廣義二階速度滑移邊界條件，并基于Succi的邊界處理方法，推導(dǎo)出壁面滑移速度和反彈比例系數(shù)的計(jì)算公式，建立了適用于過渡區(qū)微尺度氣體流動(dòng)的格子Boltzmann模型；以過渡區(qū)的微尺度Poiseuille流動(dòng)為例進(jìn)行數(shù)值分析，以驗(yàn)證本文建立的計(jì)算模型的正確性.1過渡區(qū)微尺度氣體流動(dòng)的格子Boltzmann模型1.1格子Boltzmann模型格子Boltzmann方程可以從描述稀薄氣體流動(dòng)的連續(xù)Boltzmann方程得到，其一般形式為

fα（x+eαδt，t+δt）-fα（x，t）=

-1τ[fα（x，t）-fα，eq（x，t）]+δtFα（x，t），（1）

式中： fα（x，t）表示t時(shí)刻在x處的流體粒子在α方向上的速度分布函數(shù)，可簡記為fα；

fα，eq（x，t）表示t時(shí)刻在x處的流體粒子在α方向上的平衡態(tài)分布函數(shù)，可簡記為fα，eq；

Fα（x，t）表示t時(shí)刻在x處的流體粒子在α方向上的外力，可簡記為Fα；

τ表示無量綱的松弛時(shí)間；

δt表示時(shí)間步長.

離散速度模型采用文獻(xiàn)[19]中提出的DdQm系列模型，其中， d表示空間維數(shù)， m表示離散速度個(gè)數(shù).

對(duì)于二維問題，一般采用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）

式中：

eα表示流體粒子在α方向上的離散速度；

c=δx/δt，

其中：

δx和δt分別表示網(wǎng)格步長和時(shí)間步長.

相應(yīng)的平衡態(tài)分布函數(shù)為

fα，eq=ρωα1+eα，uc2s+（eαu）22c4s-u22c4s，（3）

式中：

cs=1/3表示格子聲速；

ωα表示權(quán)系數(shù)， α=0，1，2，…，8；

ω0=4/9；

ωi=1/9（i=1，2，3，4），

ωi=1/36（i=5，6，7，8）；

ρ表示宏觀的流體密度；

u表示宏觀的流體速度.

外力項(xiàng)Fα（x，t）的模型[20]為

Fα=ωαρ1-12τeαuc2s+（eαu）eαc4sa，（4）

式中：

a表示外力加速度.

宏觀的流體密度和流體速度可通過下式計(jì)算，

ρ=∑αfα，

u=1ρ∑αfαeα+12δta.（5）1.2無量綱松弛時(shí)間與Knudsen數(shù)的關(guān)系從動(dòng)理學(xué)理論可知，氣體分子平均自由程為

λ=2μ/（ρ），

其中：

μ表示動(dòng)力粘性系數(shù)；

=8RT/π表示氣體分子平均速度；

R表示氣體常數(shù)；

T表示溫度.

另一方面，氣體分子的平均自由程還可以表示為λ=τ.

在格子Boltzmann模型中，格子速度表示為

c=χRT，

其中： χ是與模型相關(guān)的常數(shù)，對(duì)于D2Q9模型，

χ=3.

此外，需要注意的是Boltzmann方程離散時(shí)，需要對(duì)τ進(jìn)行修正，修正量為-0.5，由此可得無量綱松弛時(shí)間與Knudsen數(shù)（Kn）的關(guān)系為

τ=12+6π4 KnΔ，（6）

式中：

Δ表示計(jì)算網(wǎng)絡(luò)， Δ=1/L，其中， L為流動(dòng)特征長度.

對(duì)于微尺度氣體流動(dòng)，氣體在流經(jīng)固體壁面時(shí)，在壁面附近會(huì)產(chǎn)生一個(gè)厚度為分子平均自由程量級(jí)的Knudsen邊界層，在該區(qū)域內(nèi)分子間的碰撞頻率大大降低，分子運(yùn)動(dòng)遠(yuǎn)遠(yuǎn)沒有達(dá)到局部熱力學(xué)平衡狀態(tài)，應(yīng)力與應(yīng)變率不再滿足通常的線性關(guān)系.為使格子Boltzmann模型能夠描述Knudsen層內(nèi)的非線性速度分布，通過引入Knudsen層的有效粘性系數(shù)μe對(duì)模型進(jìn)行修正，在Knudsen層內(nèi)， μe可表示為

μe=μΨ，（7）

式中：

Ψ表示修正函數(shù).

Beskok和Karniadakis建議在過渡區(qū)采用依賴于Knudsen數(shù)的修正函數(shù)，其表達(dá)式為[21]

Ψ（Kn）=11+rKn，（8）

式中：

r表示稀薄參數(shù).

Vasilis[22]采用DSMC方法，對(duì)過渡區(qū)的Couette流和Poiseuille流進(jìn)行數(shù)值模擬，指出當(dāng)r=1.5時(shí)，此修正函數(shù)的計(jì)算結(jié)果在過渡區(qū)與DSMC的計(jì)算結(jié)果吻合很好.

Knudsen層對(duì)稀薄氣體流動(dòng)的影響可通過修正函數(shù)對(duì)粘性系數(shù)的修正，將式（6）引入到格子Boltzmann模型中，以有效粘性系數(shù)μe代替粘性系數(shù)μ，則可得考慮Knudsen層效應(yīng)后，無量綱松弛時(shí)間與Knudsen數(shù)的關(guān)系為

τ=12+6π4 KnΨ（Kn）Δ.（9）1.3微尺度氣體流動(dòng)的邊界條件對(duì)于微尺度氣體流動(dòng)，氣體在邊界上存在速度滑移現(xiàn)象，必須采用能反映速度滑移效應(yīng)的邊界格式.本文采用Succi提出的混合邊界條件[23]，即將無滑移的反彈格式與無窮滑移的鏡面反彈格式組合起來的混合格式，簡記為BSR格式.

考慮如圖1所示的氣體在二維通道中的流動(dòng)，假設(shè)密度ρ為常數(shù)， y方向的速度分量為0，并且/x=0（表示任一物理量），外力加速度a與壁面平行，采用半步長反彈的BSR格式，即壁面置于j=1/2及j=n+1/2處.

圖1二維通道中D2Q9模型的邊界布置

Fig.1Boundary setting of D2Q9 model in

two dimensional channel

在時(shí)刻t，首先執(zhí)行碰撞過程

f ′α（x，t）=fα（x，t）+

1τ[fα（x，t）-fα，eq（x，t）]+δtFα，（10）

式中：

f′α（x，t）表示碰撞后的速度分布函數(shù).

對(duì)于1

fα（x+eαδt，t+δt）=f ′α（x，t）.（11）

對(duì)于j=1的格點(diǎn)，只有f0、f1、f3、f4、f7和f8可以通過遷移過程得到， f2、f5和f6需要根據(jù)邊界條件確定；

對(duì)于j=n的格點(diǎn)， f0、f1、f2、f3、f5和f6可以通過遷移過程得到， f4、f7和f8需要根據(jù)邊界條件確定，采用BSR格式，其表達(dá)式分別為

f2=f′4，

f5=rbf′7+（1-rb）f′8，

f6=rbf′8+（1-rb）f′7，

f4=f′2，

f7=rbf′5+（1-rb）f′6，

f8=rbf′6+（1-rb）f′5，（12）

式中：

rb∈[0，1]表示反彈比例系數(shù)，即粒子在與壁面作用時(shí)沿原路彈回所占的比例.

對(duì)于定常外力驅(qū)動(dòng)下的Poiseuille流動(dòng)， BSR格式在邊界上產(chǎn)生一個(gè)無量綱滑移速度（采用中心線速度進(jìn)行無量綱化）[24]

Us=2（1-rb）rb（2τ-1）Δ+

43（2τ-1）2Δ2-Δ2，（13）

式中：

Us表示由邊界條件確定的滑移速度，

Us=us/uc，

其中， uc表示中心線上的最大速度.

根據(jù)式（9）、（13），Us還可采用Knudsen數(shù)表示為

Us=1-rbrb6πKnΨ（Kn）+

2πKn2Ψ2（Kn）-Δ2.（14）

考慮的二階滑移邊界條件為

us=C1λunwall-C2λ2unwall，（15）

式中：

C1、C2表示依賴氣體與壁面相互作用的參數(shù).

Knudsen層對(duì)二階滑移邊界條件的影響可通過修正函數(shù)對(duì)粘性系數(shù)的修正，由式（7）引入到微尺度氣體流動(dòng)的邊界條件模型中，廣義二階滑移邊界條件為

us=C1λΨ（Kn）unwall-C2λ2Ψ2（Kn）unwall.（16）

利用廣義二階滑移邊界條件，可得Poiseuille流動(dòng)在壁面上的無量綱滑移速度（采用中心線速度進(jìn)行無量綱化）為

Us=4C1KnΨ（Kn）+8C2Kn2Ψ2（Kn）.（17）

為實(shí)現(xiàn)上述滑移速度，由式（14）、（17）可知， BSR格式中的參數(shù)rb必須選擇為

rb=1+16πΔ2KnΨ（Kn）+4C1+

（8C2-2π）KnΨ（Kn）-1，（18）

式中：

C1=0.818 3；

C2=0.653 1.2數(shù)值計(jì)算管道的稀薄氣體流動(dòng)廣泛應(yīng)用于各種技術(shù)，如真空管道、微機(jī)電系統(tǒng)、空天飛行器等.在管道兩端以外力驅(qū)動(dòng)的Poiseuille流動(dòng)，引起了許多科學(xué)家和技術(shù)人員的研究興趣.由于對(duì)這類流動(dòng)的研究比較透徹，所以理論及試驗(yàn)研究結(jié)果可以作為檢驗(yàn)新模型和新算法的參考基準(zhǔn).兩個(gè)平行平板之間的微尺度Poiseuille流動(dòng)如圖1所示. Karniadakis采用局部平均速度對(duì)速度分布進(jìn)行無量綱化，得到了僅依賴于Knudsen數(shù)和y的無量綱速度剖面的表達(dá)式[25]為

U（y，Kn）=-y2L2+yL+Kn1-bKn16+Kn1-bKn，（19）

式中：

y表示距通道下壁面的距離.

當(dāng)b=-1時(shí)，用式（19）可導(dǎo)出較為準(zhǔn)確的適用于過渡區(qū)流動(dòng)的速度分布.

本文計(jì)算的Knudesn數(shù)范圍為： 0.1≤Kn≤10，流動(dòng)處于過渡區(qū).計(jì)算時(shí)取a=10-3， Δ=1/100.入口及出口采用周期性邊界條件，上下壁面采用BSR格式，稀薄參數(shù)分別取1.50和1.64.

圖2給出了不同Knudsen數(shù)下， Poiseuille流動(dòng)的無量綱速度剖面.

由圖2可知，當(dāng)稀薄參數(shù)取1.5時(shí)，在Knudsen數(shù)較小時(shí)，計(jì)算得到的無量綱速度剖面與Karniadakis給出的無量綱速度剖面吻合較好，在Knudsen數(shù)較大時(shí)，計(jì)算得到的邊界上的無量綱滑移速度大于Karniadakis給出的無量綱滑移速度，而中心線上計(jì)算得到的無量綱速度小于Karniadakis給出的無量綱速度.當(dāng)稀薄參數(shù)取 1.64時(shí)，計(jì)算得到的無量綱速度剖面在整個(gè)過渡區(qū)內(nèi)與Karniadakis給出的無量綱速度剖面吻合很好.前面的分析可知，微尺度氣體流動(dòng)會(huì)在固體邊界上產(chǎn)生Knudsen邊界層，且隨著Knudsen數(shù)的增加，Knudsen邊界層的厚度隨之增加，從而需要采用修正函數(shù)對(duì)Knudsen邊界層內(nèi)的粘性系數(shù)進(jìn)行修正，否則會(huì)高估邊界上的滑移速度.當(dāng)稀薄參數(shù)取1.50時(shí)，在Knudsen數(shù)較大時(shí)，存在修正不足的現(xiàn)象，從而導(dǎo)致邊界上的無量綱滑移速度大于Karniadakis給出的無量綱滑移速度，由于流量是固定的，從而中心線上的無量綱速度小于Karniadakis給出的無量綱速度.當(dāng)稀薄參數(shù)取1.64時(shí)，修正函數(shù)能夠?qū)nudsen邊界層內(nèi)的粘性系數(shù)進(jìn)行很好的修正，從而使得計(jì)算得到的無量綱速度剖面與Karniadakis給出的無量綱速度剖面吻合很好.

（a） Kn=0.1（b） Kn=0.4（c） Kn=0.8（d） Kn=2.0（e） Kn=6.0（f） Kn=10.0圖2微尺度Poiseuille流動(dòng)的無量綱速度剖面

Fig.2Dimensionless velocity profile of microscale Poiseuille flows由圖2還可以看出，在過渡區(qū)內(nèi)，當(dāng)Knudsen數(shù)固定時(shí)，由于壁面粘性作用的影響，壁面上的無量綱速度最小，而中心線上的無量綱速度最大，微尺度Poiseuille流動(dòng)的無量綱速度剖面基本上仍保持為拋物線形狀.隨著Knudsen數(shù)的增加，邊界上的無量綱滑移速度增大，而中心線上的無量綱速度減小.這主要是由于壁面上的無量綱滑移速度與Knudsen數(shù)成正比， Knudsen數(shù)越大，邊界上的無量綱滑移速度也越大.此外，隨著Knudsen數(shù)的增加， Knudsen層內(nèi)的有效粘性系數(shù)減小，同時(shí)由于邊界上的無量綱滑移速度增大，從而導(dǎo)致了中心線上的無量綱速度減小.

為進(jìn)一步研究稀薄參數(shù)為1.64的適用性，取a=0.01和a=0.10時(shí)，計(jì)算微尺度氣體在過渡區(qū)內(nèi)的流動(dòng)，結(jié)果如圖3所示.

由圖3可以看出，當(dāng)a=0.01， 0.10、稀薄參數(shù)取1.64時(shí)，數(shù)值計(jì)算得到的微尺度Poiseuille流動(dòng)的無量綱速度剖面仍然與Karniadakis給出的無量綱速度剖面吻合很好.

（a） Kn=0.1（b） Kn=10.0圖3不同外力下的微尺度Poiseuille流動(dòng)的無量綱速度剖面

Fig.3Dimensionless velocity profile of microscale Poiseuille flows for different body forces3結(jié)束語本文基于標(biāo)準(zhǔn)格子Boltzmann模型，建立了適用于過渡區(qū)微尺度氣體流動(dòng)的格子Boltzmann模型，并對(duì)過渡區(qū)的微尺度Poiseuille流動(dòng)進(jìn)行數(shù)值模擬.通過數(shù)值計(jì)算發(fā)現(xiàn)，稀薄參數(shù)影響壁面邊界上的滑移速度，當(dāng)稀薄參數(shù)取1.50時(shí)，計(jì)算得到的壁面上的無量綱滑移速度偏大；當(dāng)稀薄參數(shù)取1.64時(shí)，計(jì)算得到的無量綱剖面與Karniadakis給出的無量綱速度剖面吻合很好，并且在不同的外力作用下仍可以獲得很好的模擬結(jié)果，從而驗(yàn)證了本文計(jì)算方法及計(jì)算程序的正確性.參考文獻(xiàn)：[1]祝兵，宋隨弟，譚長建. 三維波浪作用下大直徑圓柱繞流的數(shù)值模擬[J]. 西南交通大學(xué)學(xué)報(bào)，2012，47（2）： 224229.

ZHU Bing， SONG Suidi， TAN Changjian. Numerical simulation for diffraction around largediameter circular cylinder subjected to threedimension wave[J]. Journal of Southwest Jiaotong University， 2012， 47（2）： 224229.

[2]OHWADA T， SONE T， AOKI K. Numerical analysis of the Poiseuille and thermal transpiration flows and a rarefied gas between two parallel plates[J]. Physics of Fluids A， 1989， 1（12）： 20422049.

[3]ABDOU M A. Chapmanenskogmaximum entropy method on timedependent neutron transport equation[J]. Journal of Quantitative Spectroscopy and Radiative Transfer， 2006， 101（2）： 210225.

[4]KUN X. Numerical hydrodynamics from gas kinetic theory[D]. New York： Columbia University， 1993.

[5]ARKILIC E B， SCHMIDT M A， BREUER K S. Gaseous slip flow in long microchannels[J]. Journal of Microelectromechanical Systems， 1997， 62（2）： 167178.

[6]BESKOK A. Simulation and models for gas flows in microg eometries[D]. Princeton： Princeton University， 1996.

[7]BIRD G A. Recent advances and current challenges for DSMC[J]. Computers and Mathematics with Applications， 1998， 35（1/2）： 114.

[8]FAN J， SHEN C. Statistical simulation of lowspeed rarefied gas flows[J]. Journal of Computational Physics， 2001， 167（2）： 393412.

[9]CAI C P ， BOYD I D ， FAN J， et al. Direct simulation methods for lowspeed microchannel flows[J]. Journal of Thermophysics and Heat Transfer， 2000， 14（3）： 368378.

[10]NIE X， DOOLEN G D， CHEN S. LatticeBoltzmann simulations of fluid flows in MEMS[J]. Journal of Statistical Physics， 2002， 107（1）： 279289.

[11]LIM C Y， SHU C， NIU X D， et al. Application of lattice Boltzmann method to simulate microchannel flows[J]. Physics of Fluids， 2002， 14（7）： 22992308.

[12]SHEN Q， TIAN D B， XIE C， et al. Examination of the LBM in simulation of microchannel flow in transitional regime[J]. Microscale Thermophysical Engineering， 2004， 8（4）： 423432.

[13]ZHANG Y， QIN R， EMERSON D R. Lattice Boltzmann simulation of rarefied gas flows in microchannels[J]. Physical Review E， 2005， 71（4）： 04770210477025.

[14]LEE T， LIN C L. Rarefaction and compressibility effects of the lattice Boltzmann equation method in a gas microchannel[J]. Physical Review E， 2005， 71（4）： 046706104670610.

[15]GUO Z L， ZHAO T S， SHI Y. Physical symmetry，spatial accuracy，and relaxation timeof the lattice Boltzmann equation for microgas flows[J]. Journal of Applied Physics， 2006， 99（7）： 07490310749038.

[16]ZHANG Y H， GU X J， BARBER R W， et al. Capturing Knudsen layer phenomena usinga lattice Boltzmann model[J]. Physical Review E， 2006， 74（4）： 04670410467044.

[17]ZHANG Y H， GU X J， BARBER R W， et al. Modelling thermal flow in the transitionregime using a lattice Boltzmann approach[J]. Europhysics Letters， 2007， 77（3）： 300031300035.

[18]GUO Z L， ZHANG C G， SHI B C. Lattice Boltzmann equation with multiple effective relaxation times for gaseous microscale flows[J]. Physical Review E， 2008， 77（3）： 036707103670712.

[19]QIAN Y H， DHUMIERES D， LALLEMAND P. Lattice BGK models for NavierStokes equation[J]. Europhysis Letters， 1992， 17（6）： 479484.

[20]GUO Z L， ZHENG C G， SHI B C. Discrete lattice effects on the forcing term in the lattice Boltzmann method[J]. Physical Review E， 2002， 65（4）： 04630810463086.

[21]BESKOK A， KARNIADAKIS G E. A model for flows in channels， pipes， and ducts at micro and nano scales[J]. Microscale Thermophys Engineering， 1999， 3（1）： 4377.

[22]VASILIS K， MICHALIS A N， KALARAKIS E D， et al. Rarefaction effects on gas viscosity in the Knudsen transition regime[J]. Microfluid Nanofluid， 2010， 9（4/5）： 847853.

[23]SUCCI S. Mesoscopic modeling of slip motion at fluid solid interfaces with heterogeneouscatalysis[J]. Physieal Review Letters， 2002， 89（6）： 06450210645024.

[24]GUO Z L， SHI B C， ZHAO T S， et al. Discrete effects on boundary conditions for the lattice Boltzmann equation in simulating mircoscale gas flows[J]. Physical Review E， 2007， 76（5）： 05670410567045.

[25]KARNIADAKIS G E， BESKOK A. Micro flows fundamental and simulation[M]. New York： SpringerVerlag， 2002： 90110.