微尺度下氣體在過渡區(qū)內(nèi)流動的格子Boltzmann模擬

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

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

中圖分類號： 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 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

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

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

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

本文在標準格子Boltzmann模型的基礎(chǔ)上，從氣體動理學理論和Knudsen效應出發(fā)，推導出適用于過渡區(qū)微尺度氣體流動的無量綱松弛時間與Knudsen數(shù)的關(guān)系，進而將二階速度滑移邊界條件推廣為廣義二階速度滑移邊界條件，并基于Succi的邊界處理方法，推導出壁面滑移速度和反彈比例系數(shù)的計算公式，建立了適用于過渡區(qū)微尺度氣體流動的格子Boltzmann模型；以過渡區(qū)的微尺度Poiseuille流動為例進行數(shù)值分析，以驗證本文建立的計算模型的正確性.1過渡區(qū)微尺度氣體流動的格子Boltzmann模型1.1格子Boltzmann模型格子Boltzmann方程可以從描述稀薄氣體流動的連續(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時刻在x處的流體粒子在α方向上的速度分布函數(shù)，可簡記為fα；

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

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

τ表示無量綱的松弛時間；

δt表示時間步長.

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

對于二維問題，一般采用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)格步長和時間步長.

相應的平衡態(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表示宏觀的流體速度.

外力項Fα（x，t）的模型[20]為

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

式中：

a表示外力加速度.

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

ρ=∑αfα，

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

λ=2μ/（ρ），

其中：

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

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

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

T表示溫度.

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

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

c=χRT，

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

χ=3.

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

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

式中：

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

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

μe=μΨ，（7）

式中：

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

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

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

式中：

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

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

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

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

考慮如圖1所示的氣體在二維通道中的流動，假設(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

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

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

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

式中：

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

對于1

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

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

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

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ù)，即粒子在與壁面作用時沿原路彈回所占的比例.

對于定常外力驅(qū)動下的Poiseuille流動， BSR格式在邊界上產(chǎ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層對二階滑移邊界條件的影響可通過修正函數(shù)對粘性系數(shù)的修正，由式（7）引入到微尺度氣體流動的邊界條件模型中，廣義二階滑移邊界條件為

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

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

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

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

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

式中：

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

當b=-1時，用式（19）可導出較為準確的適用于過渡區(qū)流動的速度分布.

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

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

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

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

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

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

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

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