LBM中基于半程反彈的統(tǒng)一邊界條件研究

凌風(fēng)如，張超英，陳燕雁，覃章榮

(廣西師范大學(xué) 廣西多源信息挖掘與安全重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004)

近年來(lái)，格子Boltzmann方法(lattice Boltzmann method，簡(jiǎn)稱LBM)在流體力學(xué)中備受關(guān)注，它用簡(jiǎn)單的微觀動(dòng)力學(xué)模型模擬復(fù)雜的宏觀現(xiàn)象，具有算法簡(jiǎn)單、邊界實(shí)現(xiàn)容易、并行度高等優(yōu)點(diǎn)[1-6]。LBM現(xiàn)已廣泛應(yīng)用于流體力學(xué)的眾多方面，如湍流、多相流、空氣動(dòng)力學(xué)、磁流體力學(xué)以及一些化學(xué)反應(yīng)和具有復(fù)雜邊界的流體流動(dòng)。對(duì)于計(jì)算流體力學(xué)中的問(wèn)題，邊界處理極其重要，它將影響模擬結(jié)果的精度和穩(wěn)定性，因此邊界處理一直是人們研究的熱點(diǎn)。

反彈格式是處理靜止無(wú)滑移平直邊界的一類常用且簡(jiǎn)單的方法。在該格式中，邊界點(diǎn)分布函數(shù)的流出方向被簡(jiǎn)單地設(shè)定為流入方向的反向[7-8]。在LBM模擬中，存在2種類型的反彈格式，即全程反彈和半程反彈[9]。在全程反彈方案中，物理邊界位于格子點(diǎn)上；在半程反彈方案中，物理邊界位于格子點(diǎn)中間，前者只有一階精度，后者在空間上具有二階精度[7]。但是，反彈格式的處理方法并不能精準(zhǔn)地模擬曲線邊界[10]。為了提高基于反彈方案的階梯逼近格式的計(jì)算精度，人們提出曲線邊界處理方法以解決復(fù)雜的幾何邊界問(wèn)題。第一種曲線邊界處理方法是由Filippova和H?nel提出的(以下簡(jiǎn)稱為FH格式)[11-12]。該方法采用一種虛平衡分布函數(shù)，基于分布函數(shù)插值構(gòu)建了具有二階精度的邊界條件。然而，它的數(shù)值穩(wěn)定性較差，特別是當(dāng)近邊界流體點(diǎn)非常接近固壁邊界時(shí)。為了提高FH格式的數(shù)值穩(wěn)定性，Mei等[13-14]提出了一種改進(jìn)的曲線邊界處理方法。Bouzidi等[15]提出了一種具有二階精度的曲線邊界處理方法，該方法結(jié)合了流體點(diǎn)分布函數(shù)的線性或二次插值和反彈技術(shù)。Lallemand等[16]改進(jìn)了這種方法，并將其擴(kuò)展到運(yùn)動(dòng)邊界中。Yu等[17]通過(guò)使用2次順序插值的方法提出一種統(tǒng)一的曲線邊界處理方法，解決了以上邊界處理的不連續(xù)性問(wèn)題。Tao等[18]提出一種基于單點(diǎn)的曲線邊界處理方法，該方法結(jié)合了插值和非平衡態(tài)外推技術(shù)，能充分實(shí)現(xiàn)LBM的并行計(jì)算優(yōu)勢(shì)。后來(lái)，學(xué)者們又進(jìn)一步改進(jìn)了這些邊界處理方法[19-21]。

在LBM模擬中，使用基于插值或外推的曲線邊界處理方法時(shí)，系統(tǒng)的質(zhì)量守恒通常難以保證，在模擬的每個(gè)時(shí)步都存在質(zhì)量的損失或增加，被稱為“質(zhì)量泄漏”。Lallemand等[16]認(rèn)為使用插值的方法會(huì)打破曲線邊界附近的質(zhì)量守恒。為了解決質(zhì)量泄漏的問(wèn)題，學(xué)者們提出基于流量的有限體積格式的曲線邊界處理方法，這無(wú)疑增加了LBM計(jì)算的復(fù)雜性[22]。隨后，學(xué)者們提出了多種質(zhì)量守恒的曲線邊界處理方法以達(dá)到避免質(zhì)量泄漏和提高數(shù)值穩(wěn)定性的效果，但這些邊界處理方法仍存在輕微的質(zhì)量泄漏問(wèn)題[23]。不過(guò)，在實(shí)際的LBM模擬中考慮到計(jì)算精度、易操作性和可靠性這3個(gè)因素，在求解曲線邊界問(wèn)題上，插值技術(shù)仍受到大多數(shù)學(xué)者們的青睞。基于此，本文提出一種質(zhì)量守恒的統(tǒng)一曲線邊界處理方法，并使用經(jīng)典算例分別驗(yàn)證該方法的計(jì)算精度、數(shù)值穩(wěn)定性和質(zhì)量守恒。

本文內(nèi)容安排如下：第1節(jié)簡(jiǎn)要回顧單弛豫的Bhatnagar-Gross-Krook (BGK)近似的格子Boltzmann方法的基本原理和文獻(xiàn)中常用的曲線邊界處理方法；第2節(jié)介紹本文提出的一種質(zhì)量守恒的統(tǒng)一曲線邊界處理新方法；第3節(jié)以Poiseuille流為例，驗(yàn)證本文提出的統(tǒng)一曲線邊界條件的有效性，并與常用的邊界條件進(jìn)行比較；第4節(jié)總結(jié)本文提出的曲線邊界方法的特點(diǎn)。

1 格子Boltzmann方法和曲線邊界條件

單弛豫時(shí)間近似的LBM模型(LBGK)的演化方程為[7]：

(1)

(2)

(3)

標(biāo)準(zhǔn)的LBM演化過(guò)程分為碰撞和流動(dòng)2個(gè)步驟，分別是：

(4)

(5)

(6)

在圖1中，當(dāng)q取值為0.5或1時(shí)，粒子分布函數(shù)以速度ei從近邊界流體點(diǎn)xf處流動(dòng)到物理邊界，接著以速度e-i反彈到xf處，這就是標(biāo)準(zhǔn)的反彈格式；但當(dāng)q取其他值時(shí)，離開(kāi)近邊界流體點(diǎn)xf處的粒子經(jīng)歷反彈后，并未回到xf處，因此需要專門的處理。

圖1 規(guī)則格子和曲線邊界的一維投影Fig.1 1D projection of regular lattice and curved boundary

(7)

(8)

其中，uf=u(xf，t)是近邊界流體點(diǎn)xf處的流體速度，usf是一個(gè)待定的虛擬速度，c=δx/δt=1。Filippova和H?nel考慮到緩慢流動(dòng)因素，對(duì)uf在xw處作Taylor展開(kāi)后進(jìn)行化簡(jiǎn)，提出用以下關(guān)系式來(lái)確定usf和χ：

(9)

(10)

(11)

(12)

當(dāng)弛豫時(shí)間τ取值較小時(shí)，Mei等改進(jìn)的公式可獲得較好的穩(wěn)定性，但是當(dāng)τ>1.25時(shí)其仍然可能存在發(fā)散的問(wèn)題[24]。除此之外，Mei等[25]還將改進(jìn)的格式(以下簡(jiǎn)稱為MLS格式)推廣到三維空間應(yīng)用。

對(duì)于無(wú)滑移的物理邊界，利用反彈的方法進(jìn)行處理是最簡(jiǎn)單和容易實(shí)現(xiàn)的，因此在LB模擬中被廣泛使用。Bouzidi等[15]提出了一種結(jié)合反彈和插值的曲線邊界處理方法(以下簡(jiǎn)稱BFL格式)。該格式的計(jì)算公式為：

(13)

(14)

在之后的研究中，Lallemand和Luo[16]應(yīng)用相同的方法改進(jìn)了BFL格式，即在引用插值公式時(shí)，把所有的分布函數(shù)統(tǒng)一到同一個(gè)時(shí)間層次上，同時(shí)考慮碰撞和流動(dòng)步驟(以下簡(jiǎn)稱LL格式)。LL格式中邊界未知分布函數(shù)的計(jì)算公式為：

(15)

(16)

(17)

可以證明上述線性插值得到的YMS1格式具有二階精度。同時(shí)可使用高階插值得到更高精度的YMS2格式[24]：

(18)

從式(17)和(18)可以看出，YMS1格式和YMS2格式在處理曲線邊界時(shí)不需要區(qū)分q的大小。然而，與其他格式相比，YMS格式需要涉及更多的操作和臨近邊界流體點(diǎn)的分布函數(shù)，這在某種程度上會(huì)影響LBM數(shù)值模擬的并行計(jì)算優(yōu)勢(shì)。

Tao等[18]提出了一種基于單點(diǎn)的二階精度的曲線邊界處理方法(以下簡(jiǎn)稱TAO格式)。TAO格式邊界處未知分布函數(shù)的計(jì)算公式為：

(19)

從其計(jì)算公式可以看出，該格式僅利用近邊界流體點(diǎn)xf處的分布函數(shù)。TAO格式使用了插值、非平衡外推和時(shí)空近似等方法，得到了一種基于單點(diǎn)的二階精度的曲線邊界處理方法，該方法能夠充分發(fā)揮LBM數(shù)值模擬的并行計(jì)算優(yōu)勢(shì)。

2 基于半程反彈的統(tǒng)一邊界條件

(20)

最后利用t+δt時(shí)刻的分布函數(shù)f-i(xv，t+δt)和f-i(x′f，t+δt)插值，即可求得邊界處待求的分布函數(shù)f-i(xf，t+δt)。

(21)

代入式(20)得f-i(xv，t+δt)，再進(jìn)行插值得：

(22)

即可求得：

(23)

(24)

(25)

其中，α是一個(gè)與τ有關(guān)的修正因子，用于減小插值帶來(lái)的誤差，由擬合Poiseuille流的數(shù)據(jù)確定，形式如下：

(26)

圖2 HBI格式示意圖Fig.2 Schematic of HBI model

3 模擬計(jì)算結(jié)果

在本文的研究中，通過(guò)2種經(jīng)典的Poiseuille流算例，一種是常規(guī)的Poiseuille流，另一種是給流場(chǎng)施加恒定重力，分別驗(yàn)證本文提出的HBI格式的計(jì)算精度、數(shù)值穩(wěn)定性和質(zhì)量守恒。然后將使用HBI格式獲得的結(jié)果與FH格式、MLS格式、BFL格式、LL格式、YMS1格式、YMS2格式和TAO格式分別進(jìn)行比較。

3.1 驗(yàn)證邊界條件的精度和穩(wěn)定性

為了驗(yàn)證HBI格式計(jì)算精度和數(shù)值穩(wěn)定性，我們采用FH格式、MLS格式、BFL格式、LL格式、YMS1格式、YMS2格式、TAO格式及HBI格式處理不同位置的固壁邊界對(duì)Poiseuille流進(jìn)行模擬。流場(chǎng)左右使用周期邊界，以恒定壓力作為體力驅(qū)動(dòng)流場(chǎng)內(nèi)分布函數(shù)運(yùn)動(dòng)，其實(shí)現(xiàn)形式如下：

(27)

(28)

其中，ν為粘滯系數(shù)，h為流場(chǎng)高度，b=y/h。由圖3可得，使用HBI格式模擬的結(jié)果u(y)與理論值uex(y)精確重合，證明本文提出的統(tǒng)一曲線邊界條件是可靠的。

圖3 使用HBI格式模擬Poiseuille流的結(jié)果Fig.3 Simulations of Poiseuille flow using HBI model

為了定量地對(duì)比HBI和其他曲線邊界條件，定義u(y)的相對(duì)L2范數(shù)誤差(L2-norm error)：

(29)

圖4 對(duì)壓力驅(qū)動(dòng)的穩(wěn)定管道流使用不同的邊界處理方法，E2隨q變化情況Fig.4 E2 as a function of q using different boundary treatments for steady state pressure-driven channel flow

3.2 驗(yàn)證邊界條件的質(zhì)量守恒

通常，在LBM模擬中，由于使用插值邊界條件，系統(tǒng)的質(zhì)量難以做到完全守恒。在模擬的過(guò)程中，每個(gè)時(shí)步都存在質(zhì)量泄漏的問(wèn)題。在許多情況下，由于質(zhì)量泄漏非常小并且最終能夠在足夠大的時(shí)間步長(zhǎng)內(nèi)接近零，因此人們往往容易忽略其對(duì)模擬結(jié)果的影響。但是，在某些情況下，這種質(zhì)量泄漏會(huì)持續(xù)增加，從而嚴(yán)重影響模擬結(jié)果的準(zhǔn)確性。

包含重力的Poiseuille流就是一個(gè)簡(jiǎn)單但非常有用的例子[26]。設(shè)流場(chǎng)的長(zhǎng)和寬分別為100和50格子單位，重力為G=-ρgj，g=0.001，q分別取值0.2和0.8，流場(chǎng)左右使用周期邊界條件[27]，計(jì)算結(jié)果如圖5所示，利用FH格式、MLS格式、BFL格式、LL格式、YMS1格式、YMS2格式和TAO格式的曲線邊界處理方法，系統(tǒng)的總質(zhì)量會(huì)以恒定速率持續(xù)下降，而使用本文提出的HBI格式計(jì)算的總質(zhì)量基本保持不變，實(shí)現(xiàn)了系統(tǒng)的質(zhì)量守恒。

圖5 使用不同的曲線邊界處理方法，重力驅(qū)動(dòng)下系統(tǒng)質(zhì)量隨時(shí)間變化情況Fig.5 System mass changes with time using different boundary treatments for steady gravity-driven flow

4 結(jié)語(yǔ)

本文提出了一種基于半程反彈的統(tǒng)一曲線邊界處理新方法，并將該方法用于對(duì)2種不同的Poiseuille流進(jìn)行數(shù)值模擬，以系統(tǒng)地測(cè)試其計(jì)算精度和數(shù)值穩(wěn)定性，檢驗(yàn)其是否存在質(zhì)量泄漏問(wèn)題。數(shù)值模擬結(jié)果表明，與常用的曲線邊界處理方法相比，本文提出的邊界處理方法具有以下優(yōu)點(diǎn)：①計(jì)算公式統(tǒng)一，保證了邊界處理的連續(xù)性；②增強(qiáng)了數(shù)值計(jì)算的穩(wěn)定性并提高了計(jì)算精度；③應(yīng)用半程反彈原理，在q取值0.5時(shí)，能夠恢復(fù)半程反彈格式；④滿足曲線邊界上的質(zhì)量守恒。本文提出的邊界處理方法適用于無(wú)滑移邊界，接下來(lái)進(jìn)一步研究將其推廣到運(yùn)動(dòng)曲線邊界中。