顧娟 黃榮宗 劉振宇 吳慧英

(上海交通大學機械與動力工程學院,上海 200240)

一種滑移區(qū)氣體流動的格子Boltzmann曲邊界處理新格式?

顧娟 黃榮宗 劉振宇 吳慧英?

(上海交通大學機械與動力工程學院,上海 200240)

(2016年12月18日收到;2017年4月5日收到修改稿)

針對滑移區(qū)復雜氣-固邊界存在速度滑移現(xiàn)象,提出了一種基于格子Boltzmann方法的非平衡態(tài)外推與有限差分相結合的曲邊界處理新格式.該格式具有可考慮實際物理邊界與網(wǎng)格線偏移量的優(yōu)勢,較傳統(tǒng)half-way DBB(di ff usive bounce-back)格式更能準確反映實際邊界情況,同時還可獲取壁面處氣體宏觀量及其法向梯度等信息.采用本文所提曲邊界處理格式模擬分析了滑移區(qū)氣體平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉Couette流問題.研究結果表明,采用曲邊界處理新格式所得結果與理論值以及文獻結果符合良好,適用于滑移區(qū)氣體流動的復雜邊界處理,且比half-way DBB格式具有更高的精度,較修正DBB格式具有更好的適應性.

曲邊界,格子Boltzmann方法,滑移區(qū),微氣體流動

1 引言

隨著微機電系統(tǒng)(micro-electro-mechanical system,MEMS)技術在生物醫(yī)療、信息通訊及國防軍工等諸多領域的逐步應用,微尺度條件下新的物理現(xiàn)象引起了越來越多學者的重視[1,2].其中對于微尺度氣體流動[3,4]問題,因氣體分子平均自由程λ與特征尺寸H兩者數(shù)量級相當,氣體稀薄效應已不可忽略.氣體稀薄程度可用Knudsen數(shù)表示,Knudsen數(shù)(Kn)的定義為Kn=λ/H.滑移區(qū)(0.001 6 Kn 6 0.1)氣體在氣-固邊界處將出現(xiàn)明顯的速度滑移現(xiàn)象,基于連續(xù)介質假定的Navier-Stokes方程需結合滑移邊界條件方可描述該問題.有別于傳統(tǒng)計算流體動力學方法,格子Boltzmann(LB)方法不受連續(xù)介質假定限制,具有微觀粒子特性和介觀物理背景,成為研究微尺度氣體流動問題的理想數(shù)值方法[5?7].

目前,針對常規(guī)尺度曲邊界處理的研究已有較多文獻報道[8?22],其中非平衡外推格式[8]已被證明可有效處理復雜曲邊界問題[13,20],前人也對Ladd[10]提出的考慮邊界速度的half-way反彈格式進行了修正以推廣到曲面邊界條件.Yin和Zhang[15]基于Ladd[10]的half-way反彈格式將虛擬節(jié)點(實際物理邊界附近流體格點與固體格點連線的中點)處速度取代原格式中的壁面速度,提高了模擬復雜邊界問題的精度.在此基礎上,史冬巖等[16]對離流體格點較近的虛擬節(jié)點速度進行修正,進一步提高了計算精度.盡管以上方法在常規(guī)尺度曲邊界處理中表現(xiàn)出較好的優(yōu)越性,但由于其不能反映因微尺度氣體-壁面相互作用引起的邊界速度滑移,無法適用于微尺度滑移區(qū)流動問題中復雜曲邊界的處理.Szalmas[23]采用貼體曲坐標下SRB(specular-re fl ection-bounceback)格式有效地模擬了微尺度復雜氣-固邊界速度滑移問題,但貼體網(wǎng)格需要求解Jacobi變換矩陣,增加了編程復雜度.另一類處理曲邊界的方法是采用笛卡爾網(wǎng)格系統(tǒng),不需要進行網(wǎng)格變換,方便處理各種復雜物理邊界,但問題在于網(wǎng)格邊界與物理邊界不重合,因此邊界條件的處理尤為重要[24].針對這個問題,Guo等[25]采用half-way DBB格式處理同心微圓柱面旋轉Couette流,但該方法僅當虛擬節(jié)點與固體壁面重合時,才能體現(xiàn)其計算精度優(yōu)勢.Suga[26]將插值格式分別引入到BB(bounce-back)格式和DS(di ff usescattering)格式中組合為IDBB(interpolation di ff usive bounce-back)格式以處理微尺度流動問題,但其組合系數(shù)的選取仍參考傳統(tǒng)half-way DBB格式.Tao和Guo[24]在DBB格式中引入實際物理邊界與網(wǎng)格線偏移量這一參數(shù),并應用其對模型松弛時間τq和參數(shù)r進行修正,將適用于平直邊界的DBB格式推廣至滑移曲邊界情況.當偏移量為半個格子時,該格式則還原為half-way DBB格式,但為保證參數(shù)τq和r取值的合理性,增加了滑移邊界條件中氣體-壁面相互作用參數(shù)取值的難度.Yang等[27]提出一種考慮速度滑移的浸入邊界格子Boltzmann方法,并用其模擬微氣體擠壓膜阻尼效應引起的空腔內(nèi)剛性物體振動,但僅驗證了該方法在不考慮稀薄效應情況下的準確度.目前針對滑移區(qū)曲邊界的速度滑移問題,精度高、適用性廣的方法鮮有研究報道.

本文提出了一種綜合Guo等[8]非平衡態(tài)外推及Le等[17]有限差分格式的曲邊界處理格式,該格式考慮了實際物理邊界與網(wǎng)格線的偏移量,較halfway DBB格式更能準確地描述實際物理邊界的情況.由于不涉及參數(shù)取值合理性問題,本文所提格式較修正DBB格式更能適應氣體-壁面相互作用參數(shù)的選取,同時還可獲得壁面處氣體宏觀量及其法向梯度等信息.基于本文所提格式對滑移區(qū)氣體平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉Couette流問題進行LB數(shù)值模擬研究,將模擬結果與解析解及已有文獻結果進行對比討論,以分析本文所提格式的準確度與適用性.

2 LB模型與曲邊界格式

2.1 微流動的LB基本模型

LB模型通過格點上微觀粒子分布函數(shù)的演化過程來描述流體的流動.本文采用單松弛(LBGK)模型[28],其演化方程為

其中,ci為粒子在i方向的速度,δt為時間步長,fi(x,t)為t時刻x位置處以速度ci運動的粒子的分布函數(shù),τ為無量綱松弛時間.(x,t)為i方向粒子的平衡態(tài)分布函數(shù),

Fi(x,t)為離散外力項[29]

式中,ωi為權系數(shù),cs為格子聲速,G為外力.宏觀物理量密度ρ和速度u的計算式為

對于二維問題,本文選用二維九速(D2Q9)模型,相應的離散速度為

其中,粒子遷移速率c=δx/δt(δx為格子步長),格子聲速權系數(shù)ω0=4/9,ω1=ω2=ω3=ω4=1/9,ω5=ω6=ω7=ω8=1/36.通過Chapman-Enskog分析可知,模型對應的流體壓力P=,流體的運動黏度ν=(τ?0.5)δt.

對于微尺度氣體流動問題,需在模型中引入Knudsen數(shù)的影響.由動理學理論可知,氣體分子平均自由程λ與動力黏度μ、壓力P、溫度T的關系式為結合Knudsen數(shù)的定義式Kn=λ/H及LB模型中運動黏度ν與無量綱松弛時間τ的關系式ν=c2s(τ?0.5)δt,可得模型中無量綱松弛時間與Knudsen數(shù)的關系式為

2.2微尺度流動的曲邊界處理

對于微尺度氣體流動問題,壁面處滑移邊界條件的處理至關重要.本文在Guo等[8]的非平衡態(tài)外推格式的基礎上,提出了一種復雜曲邊界上壁面滑移速度的邊界處理格式.如圖1所示的曲邊界,xf為流體格點,xs為固體內(nèi)部格點,xint為固體邊界格點.流體格點xf處i方向的未知分布函數(shù)fi(xf,t+δt)由固體格點xs處的分布函數(shù)fi(xs,t)演化得到,而fi(xs,t)采用Guo等[8]的非平衡態(tài)外推格式重構得到.為此,首先將固體格點處的分布函數(shù)分解為平衡態(tài)和非平衡態(tài)兩部分,即

式中us可由u(xint),u(xf)及u(xff)外推確定

需要指出的是,對于微尺度氣體流動問題,壁面處氣體速度ug(即u(xint))不等于固體壁面速度uw,需進一步考慮滑移速度uslip的影響.本文采用微尺度一階滑移曲邊界條件[30],

其中,uτ為氣體在壁面處的切向速度矢量,r為半徑,σv為動量協(xié)調(diào)系數(shù).對于平直邊界(r→∞)也適用,此時(11)式括號內(nèi)第二項為0,與微尺度經(jīng)典的Maxwell-Smoluchowski一階滑移直邊界條件[1]表達式一致.對于完全漫反射固體壁面,有σv=1[3?7,23?27].針對復雜曲邊界,本文采用Le等[17]提出求解壁面法向溫度梯度的有限差分格式求解壁面處流體速度的法向梯度.如圖1所示,沿壁面法向n在流體內(nèi)部選取兩個虛擬點x′和x′′,滿足|x′′?x′|=|x′?xint|=δ.為保證緊鄰x′(x′′)的四個網(wǎng)格點均為流體格點,本文取δ=1.5δx.對于虛擬點x′處的流體宏觀速度,通過其緊鄰的流體網(wǎng)格點采用雙線性插值法確定,即有

其中A1,A2,A3和A4分別為x′所在網(wǎng)格內(nèi)四塊子區(qū)域的面積(如圖1所示).采用相同方法可求得虛擬點x′′處的流體宏觀速度u′′.由虛擬點x′和x′′處的流體宏觀速度,結合(11)式,可求得壁面處氣體速度為

其中第二項即為微尺度壁面處氣體滑移速度,R為固壁的曲率半徑.

圖1 曲邊界示意圖Fig.1.Schematic illustration of the curved boundary.

3 結果與討論

為驗證本文提出的滑移區(qū)曲邊界滑移速度邊界處理格式的適用性,本節(jié)依次模擬研究平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉Couette流問題.

3.1 平直微通道Poiseuille流

微尺度Poiseuille流在MEMS中廣泛存在且有解析解,可用于驗證滑移邊界處理方法的準確度.本文模擬了網(wǎng)格線與平直微通道壁面偏移量qδx(q=0.25,0.5,0.75,1.0)不同情況下的滑移區(qū)外力驅動Poiseuille流.微通道內(nèi)氣體受到沿流動方向恒定外力Gx=10?4的驅動,進出口采用周期性邊界條件,固體邊界上的滑移速度邊界條件((11)式)采用本文提出的邊界處理格式處理.計算過程中,計算網(wǎng)格取為Nx×Ny=4×32.對于該微尺度Poiseuille流,沿截面方向速度分布的解析解為

其中H為微通道的寬度.

圖2(a)和圖2(b)分別給出了當Kn=0.0194時,采用本文格式和half-way DBB格式計算得到的截面無量綱速度分布與解析解的對比情況.由圖2可知,本文格式得到的截面無量綱速度分布和解析解在不同q情況下均符合良好;而half-way DBB格式得到的結果僅在q=0.5時才與解析解一致,這是因為half-way DBB格式基于half-way反彈格式,當偏移量為0.5δx時,half-way DBB格式中的虛擬邊界與實際物理邊界重合,固體邊界可由虛擬邊界準確刻畫,故此時計算精度最高.但當偏移量發(fā)生改變時,half-way DBB格式因默認偏移量為0.5δx,故其模擬結果并未有相應改變,仍對應q=0.5條件下的計算結果.故對于偏移量不為0.5δx,half-way DBB格式已不能再準確描述固體邊界的實際情況.

為定量比較本文格式與half-way DBB格式的計算優(yōu)勢,定義了如下誤差量:1)最大相對誤差,

2)相對誤差,

其中,(i,j)為計算域上的網(wǎng)格點,u(i,j)為流體格點上的速度值,uexact(i,j)為流體格點上的速度的解析值.表1給出了不同q條件下(Kn=0.0194),本文格式與half-way DBB格式計算結果誤差對比情況.從表中可知,如前所述half-way DBB格式僅在q=0.5時誤差才達最小值,本文格式在其余不同偏移量條件下較half-way DBB格式計算精度都有明顯提高.

圖2 (網(wǎng)刊彩色)不同q時平直微通道截面無量綱速度分布(Kn=0.0194)(a)本文格式模擬結果;(b)half-way DBB格式模擬結果Fig.2.(color online)Normalized velocity pro fi les of the aligned microchannel fl ow with di ff erent values of q(Kn=0.0194):(a)The present scheme;(b)the half-way DBB scheme.

3.2 傾斜微通道Poiseuille流

為進一步驗證本文提出的格式能有效處理滑移區(qū)復雜的幾何邊界,對圖3所示的傾斜微通道內(nèi)外力驅動Poiseuille流進行了模擬.微通道傾斜角度θ依次取為arctan(0.46),arctan(1.20)和arctan(3.72),計算網(wǎng)格取為Nx×NAB=100×100,對應的y方向格點數(shù)為Ny=NAB+Nxtanθ.對于該傾斜微通道Poiseuille流,截面速度分布解析解仍如(14)式所示,其中y代表到微通道下壁面的垂直距離,H為微通道的寬度.

圖3傾斜微通道示意圖Fig.3.Schematic illustration of the inclined microchannel.

圖4 (a)和圖4(b)分別給出了當Kn=0.0194時,采用本文格式和half-way DBB格式計算得到的截面無量綱速度分布與解析解對比情況.由圖4可知,在不同傾斜角度情況下,本文格式的模擬結果和解析解始終符合良好,而half-way DBB格式的模擬結果在微通道壁面處將顯著偏離理論值.這是因為halfway DBB格式實際物理邊界與網(wǎng)格線偏移量恒為0.5δx,使得傾斜的微通道壁面成為階梯狀.

表2給出了不同傾斜角條件下(Kn=0.0194),本文格式與half-way DBB格式計算結果誤差對比情況.從表2可知,本文格式較half-way DBB格式在模擬不同傾斜角時最大相對誤差及相對誤差均明顯降低.

圖4 (網(wǎng)刊彩色)不同傾斜角θ時,傾斜微通道截面無量綱速度分布(Kn=0.0194)(a)本文格式模擬結果;(b)half-way DBB格式模擬結果Fig.4.(color online)Normalized velocity pro fi les of the inclined microchannel fl ow with di ff erent inclined angles θ(Kn=0.0194):(a)The present scheme;(b)the half-way DBB scheme.

表2 不同傾斜角條件下本文格式與half-way DBB格式計算誤差對比(Kn=0.0194)Table 2.Computational errors of the present scheme and the half-way DBB scheme with di ff erent inclined angles θ(Kn=0.0194).

表1 不同q條件下本文格式與half-way DBB格式計算誤差對比(Kn=0.0194)Table 1.Computational errors of the present scheme and the half-way DBB scheme with di ff erent values of q(Kn=0.0194).

3.3微圓柱氣體繞流

圓柱繞流是流體力學中經(jīng)典的復雜外流問題,常被用于檢測曲邊界的處理能力.本節(jié)對微圓柱氣體繞流問題進行模擬,并將模擬結果與文獻結果[31]進行對比驗證.對于微圓柱繞流問題,Reynolds數(shù)的定義為

其中,U∞為均勻來流速度,D為微圓柱直徑.在模擬過程中,計算區(qū)域選取為40D×20D(D劃分為30個格子),微圓柱放置在離進口及上下邊界均為10D的位置.計算區(qū)域邊界條件設置如下:

圖5 (網(wǎng)刊彩色)微圓柱繞流(Re=20,Kn=0.015)(a)微圓柱附近流線;(b)微圓柱表面無量綱滑移速度分布Fig.5.(color online)Flow past a microcylinder(Re=20,Kn=0.015):(a)Streamlines around the microcylinder;(b)normalized slip velocity pro fi les on the surface of the microcylinder.

微圓柱表面則采用(11)式所示的一階滑移曲邊界條件.

圖5(a)給出了Re=20時,考慮稀薄效應(Kn=0.015)與不考慮稀薄效應(Kn=0)情況下微圓柱附近流線分布的對比.從圖中可以看出,考慮稀薄效應時回流區(qū)長度明顯縮短,這是由于氣體-壁面滑移速度使得壁面附近流體微團動能增加,逆壓梯度減小,從而分離渦變小.圖5(b)給出了Kn=0.015,Re=20時,采用本文格式及half-way DBB格式模擬得到的微圓柱表面無量綱滑移速度分布同Beskok和Karniadakis[31]采用譜元方法計算得到結果的對比.由圖可知,本文格式得到的模擬結果與譜元方法計算結果更為接近,在分離點之前幾乎完全一致,在分離點之后存在一定偏差;而half-way DBB格式得到的模擬結果整體偏小,最大滑移速度較文獻[31]的結果存在較大偏差.

3.4 同心微圓柱面Couette流

上述算例均為靜止壁面上的滑移速度邊界條件,為驗證本文提出的格式亦能處理運動壁面上的滑移速度邊界條件,本文模擬了同心微圓柱面旋轉Couette流問題.如圖6所示,內(nèi)微圓柱的半徑為R1(R1=30δx),以恒定的角速度ω逆時針旋轉;外微圓柱的半徑為R2,保持為靜止狀態(tài);內(nèi)外微圓柱的半徑比R1:R2=3:5,其表面均為完全漫反射表面(σv1=σv2=1).

對于同心微圓柱面旋轉Couette流,圓柱坐標系下的Navier-Stokes方程可簡化為[23?26,30]

其中r為半徑,uθ為切向速度.內(nèi)外圓柱壁面處的氣體速度邊界條件為

由此可得速度分布解析解為

其中系數(shù)A1,A2為

圖6同心微圓柱面Couette流示意圖Fig.6.Schematic illustration of the concentric microcylinder Couette fl ow.

圖7 給出了Kn=0.05和Kn=0.1時,采用本文提出的格式、half-way DBB格式及Tao和Guo[24]提出的修正DBB格式計算得到的徑向方向無量綱速度分布與解析解對比.由圖7可知,本文格式的模擬結果始終與解析解符合良好;而halfway DBB格式和修正DBB格式模擬結果分別在內(nèi)微圓柱附近和外微圓柱附近存在較大偏差;且Kn數(shù)越大,half-way DBB格式和修正DBB格式在邊界處的偏差越明顯.

表3給出了不同Kn數(shù)條件下,本文格式同half-way DBB格式和修正DBB格式計算結果誤差對比情況.從表中可知,本文格式在模擬不同Kn數(shù)時最大相對誤差及相對誤差最小.且由表3可知,修正DBB格式在邊界處最大相對誤差較half-way DBB格式優(yōu)勢并不明顯,但由于在其余區(qū)域符合較好,故在相對誤差方面修正DBB格式比half-way DBB格式計算精度有所提高.

此外,為檢驗新格式的收斂性能,定義了流場的計算收斂判據(jù):

其中,u(t+1)(i,j)為當前時刻流體格點的速度值,u(t)(i,j)為前一時刻流體格點的速度值.圖8給出了Kn=0.05時,本文格式同half-way DBB格式和修正DBB格式收斂過程對比情況.由圖8可知,三種邊界處理格式的計算收斂時間相近.基于上述計算精度及收斂性能分析比較,本文提出的格式能更好地處理滑移區(qū)運動曲邊界問題.

圖7 (網(wǎng)刊彩色)微圓柱間Couette流無量綱速度分布(a)Kn=0.05;(b)Kn=0.1Fig.7.(color online)Normalized velocity pro fi les of the microcylinder Couette fl ow:(a)Kn=0.05;(b)Kn=0.1.

表3 不同Kn數(shù)條件下本文格式與half-way DBB格式及修正DBB格式計算誤差對比Table 3.Computational errors of the present scheme,the half-way DBB scheme and the modi fi ed DBB scheme with di ff erent Kn numbers.

圖8 (網(wǎng)刊彩色)不同格式的收斂曲線(Kn=0.05)Fig.8.(color online)Convergence pro fi les of di ff erent schemes(Kn=0.05).

4 結論

本文基于非平衡態(tài)外推格式處理曲邊界與有限差分格式獲取曲邊界處宏觀物理量及其法向梯度相結合的思想,提出了一種適用于滑移區(qū)復雜氣-固邊界的格子Boltzmann曲邊界處理格式.該格式在處理復雜邊界時考慮了實際物理邊界與格線偏移量的影響,較half-way DBB格式具有更高的準確度;且該格式對氣體-壁面相互作用參數(shù)選取無任何限定,相較修正DBB格式有更好的適應性.基于該格式對滑移區(qū)氣體平直/傾斜微通道內(nèi)Poiseuille流、微圓柱繞流和同心微圓柱面旋轉Couette流問題進行計算分析,得到以下結論:

1)本文所提邊界處理格式能夠準確模擬網(wǎng)格線與實際物理邊界不重合的平直以及傾斜微通道外力驅動Poiseuille流;而half-way DBB格式僅在特定偏移量條件下才能準確模擬平直微通道內(nèi)Poiseuille流,且在模擬傾斜微通道內(nèi)Poiseuille流時在壁面附近較理論值存在明顯偏離;

2)通過對微圓柱氣體繞流計算分析,發(fā)現(xiàn)考慮稀薄效應時回流區(qū)長度明顯縮短;本文所提格式得到的微圓柱表面無量綱滑移速度與文獻報道結果更為接近,而half-way DBB格式得到的結果整體偏小,且最大滑移速度存在較大偏差;

3)在模擬滑移區(qū)同心微圓柱面Couette流時,本文格式預測的徑向方向無量綱速度分布與解析解符合良好,模擬精度較half-way DBB格式及修正DBB格式有明顯提高.

綜上可知,本文提出的曲邊界處理格式具有處理滑移區(qū)復雜流固邊界速度滑移問題的能力,且較half-way DBB格式和修正DBB格式分別有更高的精度和更好的適應性.此外,本文提出的曲邊界處理格式,亦可推廣處理微氣體流動傳熱問題復雜流-固邊界溫度跳躍邊界條件.

[1]Gad-el-Hak M 1999 ASME J.Fluids Eng.121 5

[2]Ho C M,Tai Y C 1998 Annu.Rev.Fluid Mech.30 579

[3]Zhang T T,Jia L,Wang Z C 2008 Chin.Phys.Lett.25 180

[4]Sun X M,He F,Ding Y T 2003 Chin.Phys.Lett.20 2199

[5]Liu C F,Ni Y S 2008 Chin.Phys.B 17 4554

[6]Tao S,Wang L,Guo Z L 2014 Acta Phys.Sin.63 214703(in Chinese)[陶實,王亮,郭照立2014物理學報63 214703]

[7]Wang Z,Liu Y,Zhang J Z 2016 Acta Phys.Sin.65 014703(in Chinese)[王佐,劉雁,張家忠2016物理學報65 014703]

[8]Guo Z L,Zheng C G,Shi B C 2002 Phys.Fluids 14 2007

[9]Mei R W,Luo L S,Shyy W 1999 J.Comput.Phys.155 307

[10]Ladd A J C 1994 J.Fluid Mech.271 285

[11]Noble D R,Torczynski J R 1998 Int.J.Mod.Phys.C 9 1189

[12]Yu Z,Yang H,Fan L S 2011 Chem.Eng.Sci.66 3441[13]Lee C H,Huang Z H,Chiew Y M 2015 Eng.Appl.Comp.Fluid 9 370

[14]Guo K K,Li L K,Xiao G,AuYeung N A,Mei R W 2015 Int.J.Heat Mass Tran.88 306

[15]Yin X W,Zhang J F 2012 J.Comput.Phys.231 4295

[16]Shi D Y,Wang Z K,Zhang A M 2014 Acta Phys.Sin.63 074703(in Chinese)[史冬巖,王志凱,張阿漫2014物理學報63 074703]

[17]Le G G,Oulaid O,Zhang J F 2015 Phys.Rev.E 91 033306

[18]Fu S C,Leung W W F,So R M C 2013 Commun.Comput.Phys.14 126

[19]Gao D Y,Chen Z Q,Chen L H,Zhang D L 2017 Int.J.Heat Mass Tran.105 673

[20]Izham M,Fukui T,Morinishi K 2011 J.Fluid Sci.Tech.6 812

[21]Nie D M,Lin J Z 2010 Chin.Phys.Lett.27 104701

[22]Fang H P,Wan R Z,Fan L W 2000 Chin.Phys.9 515

[23]Szalmas L 2007 Int.J.Mod.Phys.C 18 15

[24]Tao S,Guo Z L 2015 Phys.Rev.E 91 043305

[25]Guo Z L,Shi B C,Zheng C G 2011 Comput.Math.Appl.61 3519

[26]Suga K 2013 Fluid Dyn.Res.45 034501

[27]Yang W L,Li H X,Zhang T J,Elfadel I M 2015 ASME 2015 13th International Conference on Nanochannels,Microchannels,and Minichannels San Francisco,CA,USA,July 6–9,2015 pV001T04A001

[28]Buick J M,Greated C A 2000 Phys.Rev.E 61 5307

[29]Guo Z L,Zheng C G,Shi B C 2002 Phys.Rev.E 65 046308

[30]Sun Y H,Barber R W,Emerson D R 2005 Phys.Fluids 17 047102

[31]Beskok A,Karniadakis G E 1994 J.Thermophys.Heat Tr.8 647

PACS:47.11.Qr,47.15.Rq,02.60.CbDOI:10.7498/aps.66.114701

A new curved boundary treatment in lattice Boltzmann method for micro gas fl ow in the slip regime?

Gu JuanHuang Rong-ZongLiu Zhen-YuWu Hui-Ying?
(School of Mechanical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

18 December 2016;revised manuscript

5 April 2017)

A new curved boundary treatment in lattice Boltzmann method is developed for micro gas fl ow in the slip regime.The proposed treatment is a combination of the nonequilibrium extrapolation scheme for curved boundary with noslip velocity condition and the counter-extrapolation method for the velocity and its normal gradient on the curved boundary.Taking into consideration the e ff ect of the o ff set between the physical boundary and the closest grid line,the new treatment is proved to be more accurate than the traditional half-way di ff usive bounce-back(DBB)scheme.The present treatment is also more applicable than the modi fi ed DBB scheme because the speci fi c gas-wall interaction parameters need to be determined to ensure the validation of the modi fi ed DBB scheme.

The proposed boundary treatment is implemented to simulate the benchmark problems,which include a Poiseuille fl ow in the aligned/inclined micro-channel,a fl ow past a microcylinder and a microcylindrical Couette fl ow.The results and conclusions are summarized as follows.

1)The force-driven Poiseuille fl ow in an aligned microchannel is simulated separately with di ff erent values of wallgrid o ff set qδx(q=0.25,0.5,0.75,1.0).With the consideration of the wall-grid o ff set,the numerical results with the new boundary treatment show good agreement with the analytical solutions.However,the results obtained by using the half-way DBB scheme only accord well with the analytical solutions under the condition of a fi xed wall-grid o ff set(q=0.5).

2)To demonstrate the capability of the present treatment in dealing with gas fl ow in a more complex geometry,the force-driven Poiseuille fl ow in a micro-channel is investigated separately with di ff erent inclined angles.The present numerical results fi t well with the analytical solutions.However,the discrepancy between the results obtained with the half-way DBB scheme and the analytical solutions can be clearly observed near the inclined boundaries.

3)The gas fl ow past a microcylinder is simulated to prove that the present treatment can deal with the curved boundary.The slip velocity pro fi le along the micro cylinder periphery obtained with the present treatment accords well with the available data in the published literature.However,the results obtained with the half-way DBB scheme show lower values than the results from the published work.

4)In the simulations of the microcylindrical Couette fl ow between two coaxial rotating cylinders for di ff erent Knudsen numbers the results obtained by using the present treatment show excellent agreement with the analytical solutions.However,the results obtained with the half-way DBB scheme and the modi fi ed DBB scheme deviate obviously from the analytical solutions near the inner and outer cylindrical walls,respectively.

In summary,the new boundary treatment proposed in this work is capable of dealing with the complex gas-solid boundary in the slip regime.The new treatment has a higher accuracy than the half-way DBB scheme and shows a better applicability than the modi fi ed DBB scheme.

curved boundary,lattice Boltzmann method,slip regime,micro gas fl ow

10.7498/aps.66.114701

?國家自然科學基金(批準號:51536005,51521004)資助的課題.

?通信作者.E-mail:whysrj@sjtu.edu.cn

?2017中國物理學會Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.51536005,51521004).

?Corresponding author.E-mail:whysrj@sjtu.edu.cn