慣量松弛因子在格子Boltzmann 方法中的應(yīng)用

王穎娟, 龔光彩, 石 星, 龔子徹, 劉永超

(湖南大學(xué)土木工程學(xué)院，長(zhǎng)沙 410082)

1 引言

格子Boltzmann 方法(LBM)是近年來(lái)發(fā)展較迅速的一種新型數(shù)值計(jì)算方法，不同于傳統(tǒng)的計(jì)算流體力學(xué)方法，它是一種介于宏觀與微觀之間的介觀方法[1]，主要有以下優(yōu)點(diǎn)：算法簡(jiǎn)單、編程容易、能夠處理復(fù)雜的邊界條件、具有良好的并行性、能直接模擬有復(fù)雜幾何邊界的諸如多孔介質(zhì)等連通域流場(chǎng)，無(wú)需作計(jì)算網(wǎng)格的轉(zhuǎn)換等等.LBM 是早期的格子氣自動(dòng)機(jī)模型的拓展[2]，現(xiàn)在已經(jīng)是一種理論比較完備、模型比較成熟的數(shù)值模擬方法，目前被廣泛應(yīng)用在了多個(gè)領(lǐng)域，包括：熱效應(yīng)[3]、層流，湍流模擬[4]、復(fù)雜邊界以及動(dòng)邊界[5]等等，同時(shí)格子Boltzman 方法也可以解決傳統(tǒng)計(jì)算方法中比較難模擬的多孔介質(zhì)的問(wèn)題，將在工程材料與能源利用領(lǐng)域有著重要應(yīng)用.

慣量松弛因子一般運(yùn)用在解決N-S 方程中，在SIMPLE 算法中加入慣量松弛因子能加快程序收斂速度，提高穩(wěn)定性，并且有一定的適應(yīng)非均勻網(wǎng)格的能力[6,7].對(duì)于較為復(fù)雜的流動(dòng)，傳統(tǒng)格子Boltzmann 方法往往需要足夠多的網(wǎng)格點(diǎn)來(lái)計(jì)算流場(chǎng)，以便更好的捕捉流場(chǎng)信息，但同時(shí)也會(huì)增加計(jì)算量，使得計(jì)算效率降低.為了解決傳統(tǒng)格子Boltzmann 方法效率低的問(wèn)題，筆者首次討論了將慣量松弛因子與格子Boltzmann 方法結(jié)合，以頂蓋驅(qū)動(dòng)方腔流為校核算例，對(duì)不同慣量松弛因子下的格子Boltzmann 方法開(kāi)展數(shù)值算法研究及對(duì)應(yīng)的算例驗(yàn)證，探究慣量松弛因子在格子Boltzmann 方法中的作用及其影響.

2 慣量松弛因子與格子Boltzmann 模型

2.1 慣量松弛因子與LBGK 模型

LBGK 模型是目前應(yīng)用比較廣泛的格子Boltzmann 模型，Qian 等[8]提出的DmQn 模型(m維空間，n個(gè)離散格子速度)最為典型，為了保證各向同性、伽利略不變性及不可壓N-S 方程中速度對(duì)壓力的獨(dú)立性，本文采用LBGK 的D2Q9 模型[9]，其演化方程為

其中fk(x,t)是ck方向的粒子密度分布函數(shù)；fkeq(x,t)是在t時(shí)刻x處的平衡態(tài)分布函數(shù)；ω為松弛頻率；ck為粒子速度矢量.

其中格子速度c=δx/δt.宏觀速度、體積平均密度可以由分布函數(shù)得出

平衡態(tài)分布函數(shù)為

通過(guò)Chapman-Enskog 展開(kāi)，可以得到運(yùn)動(dòng)粘性系數(shù)與松弛頻率的關(guān)系為

LBM 方法包含兩個(gè)步驟：碰撞和遷移.

1) 碰撞步驟的迭代格式如下

在上式中引入慣量項(xiàng)Rfk，在迭代過(guò)程中，上式的迭代格式成為

其中R為慣量松弛因子.

2) 遷移步驟如下

2.2 多松弛格子Boltzmann 模型

在LBGK 模型中，當(dāng)松弛時(shí)間過(guò)大時(shí)，會(huì)出現(xiàn)數(shù)值不穩(wěn)定的情況，因此為了克服這一缺點(diǎn)，法國(guó)學(xué)者[10]提出了一個(gè)廣義格子Boltzmann 模型，即多松弛格子Boltzmann 模型(MRT-LBM).

其中f 是一個(gè)矢量，代表格子節(jié)點(diǎn)上的速度分布函數(shù)，m 代表一系列相互獨(dú)立的矩，meq是矩m 的平衡態(tài)值，M是一個(gè)正交的轉(zhuǎn)換矩陣，可把f 轉(zhuǎn)換為m

其中

松弛矩陣

3 物理模型

本文采用頂蓋驅(qū)動(dòng)流作為校核算例，它是一個(gè)經(jīng)典的不可壓縮流動(dòng)，其結(jié)果廣泛得到了驗(yàn)證.圖1 為二維方腔頂蓋驅(qū)動(dòng)流的示意圖，在實(shí)際計(jì)算中，參考文獻(xiàn)[9]中的頂蓋驅(qū)動(dòng)流的格子Boltzmann 程序，并且對(duì)采用不同慣量松弛因子的計(jì)算格式進(jìn)行代碼編譯，邊界處理采用標(biāo)準(zhǔn)反彈格式，程序收斂判據(jù)如下

其中error 為兩個(gè)相鄰時(shí)層速度的最大相對(duì)誤差；ux(i,j,t)為點(diǎn)(i,j,t)處沿x軸方向的宏觀速度；uy(i,j,t)為點(diǎn)(i,j,t)處沿y軸方向的宏觀速度；ε為一個(gè)小量.

圖1 二維頂蓋驅(qū)動(dòng)方腔流示意圖

在標(biāo)準(zhǔn)的D2Q9 格子Boltzmann 模型模擬中，頂蓋驅(qū)動(dòng)速度U=0.1，網(wǎng)格為100×100，運(yùn)動(dòng)粘度系數(shù)由雷諾數(shù)Re 定義式反算得到，各參數(shù)均為格子單位，例如對(duì)雷諾數(shù)Re 取1000，ε則取7.5E?6.筆者分別采用了LBGK 與多松弛的格式對(duì)方腔流的流函數(shù)等值線、中軸線上的速度等進(jìn)行結(jié)果對(duì)比；同時(shí)在D2Q9 的LBGK 模型中加入慣量松弛因子，并分別給出了不同慣量松弛因子下的收斂步數(shù)、流函數(shù)等值線圖、以及記錄渦心的坐標(biāo)，并將計(jì)算結(jié)果與不加慣量松弛因子的結(jié)果進(jìn)行比較分析；最后對(duì)三維的頂蓋方腔流進(jìn)行了模擬計(jì)算，研究慣量松弛因子在三維格子Boltzmann 中的作用.

4 計(jì)算結(jié)果及分析

4.1 D2Q9 的LBGK 模型與MRT-LBM 模型

筆者分別采用D2Q9 的LBGK 模型和多松弛的格子Boltzmann 模型對(duì)頂蓋驅(qū)動(dòng)方腔流進(jìn)行模擬，其中MRT-LBM 模型的程序是參考文獻(xiàn)[9]中附錄里的計(jì)算機(jī)代碼，雷諾數(shù)取1000，圖2 為兩種模型方法所求得的流函數(shù)等值線對(duì)比圖，從圖中可以看出兩種模型所求得的流函數(shù)圖趨勢(shì)基本一致，當(dāng)流體穩(wěn)定后，方腔中央、左下角和右下角都分別有一個(gè)渦.表1 為兩種模型的渦心坐標(biāo)與基準(zhǔn)解的比較，基準(zhǔn)解為文獻(xiàn)[11]中的模擬結(jié)果，可以論證兩種模型程序的正確性.

表1 數(shù)值模擬結(jié)果對(duì)比

圖3 表示兩種模型所求得的水平、垂直兩條中軸線上的速度對(duì)比圖，從圖中可以看出，兩種模型的中軸線上的速度曲線吻合的很好，基本重合一致.

圖2 LBGK 模型和MRT-LBM 模型的流函數(shù)等值線對(duì)比圖

圖3 兩種模型中軸線上的速度對(duì)比

4.2 慣量松弛因子在D2Q9 的LBGK 模型中的應(yīng)用

表2 表示雷諾數(shù)分別取400、1000 時(shí)，不加慣量松弛因子與加不同慣量松弛因子下的收斂迭代步數(shù)與計(jì)算結(jié)果偏差的對(duì)比，圖4 為收斂步數(shù)與慣量松弛因子的關(guān)系曲線圖，其中R= 0 為不加慣量松弛因子的模擬結(jié)果，從表中可以看出：當(dāng)慣量松弛因子R大于0.03 時(shí)，程序收斂步數(shù)隨著慣量松弛因子的增大在不斷的減少，相應(yīng)提高的計(jì)算效率也越來(lái)越高，最高達(dá)到了50%；且隨著雷諾數(shù)的增大，在相同的慣量松弛因子下程序所提高的效率也在增大；同時(shí)經(jīng)研究發(fā)現(xiàn)當(dāng)慣量松弛因子R取0.01 到0.02 時(shí)，我們發(fā)現(xiàn)程序的運(yùn)行結(jié)果難以達(dá)到收斂精度，但程序并不發(fā)散，這表明慣量松弛因子無(wú)論是運(yùn)用在SIMPLE 算法中還是運(yùn)用在格子Boltzmann 中都有其適宜的取值區(qū)間，因此慣量松弛因子R取0.02 時(shí)是其極限的取值下界.

表2 不同慣量松弛因子下的收斂步數(shù)與計(jì)算誤差

圖4 收斂步數(shù)與慣量松弛因子的關(guān)系曲線

圖5、圖6 分別為雷諾數(shù)取400、1000時(shí)，取不同慣量松弛因子下的水平、垂直兩條中軸線上的速度與不加慣量松弛因子的該速度進(jìn)行比較，從圖中可以看出：當(dāng)慣量松弛因子R等于0.03 時(shí)，計(jì)算結(jié)果的偏差在10%以內(nèi)；當(dāng)慣量松弛因子R等于0.05 時(shí)，計(jì)算結(jié)果的偏差在15%以內(nèi)；當(dāng)R等于0.1 時(shí)，計(jì)算結(jié)果的偏差在25%以內(nèi)；因此隨著慣量松弛因子R的增大，中軸線上的速度偏差也越來(lái)越大；而且隨著雷諾數(shù)的增大，取不同慣量松弛因子之間的誤差相對(duì)也會(huì)增大一點(diǎn)，但都在可接受范圍之內(nèi).將加了慣量松弛因子的LBGK 模型與MRT-LBM 模型進(jìn)行對(duì)比可以發(fā)現(xiàn)，當(dāng)慣量松弛因子取0.03 和0.05 時(shí)，中軸線上的速度與MRT-LBM 模型中軸線上的速度非常接近.

圖5 水平中軸線上的速度對(duì)比

圖6 垂直中軸線上的速度對(duì)比

圖7 給出了雷諾數(shù)取1000 時(shí)，不加慣量松弛因子與加不同慣量松弛因子下頂蓋驅(qū)動(dòng)流的流函數(shù)圖，從圖中可以清晰的看到：流動(dòng)穩(wěn)定后，方腔的中央都有個(gè)一級(jí)大渦，而且左下角和右下角也都分別有一個(gè)二級(jí)渦；但隨著慣量松弛因子的增大，中心渦渦心有由方腔中心逐漸向方腔中心偏右上方移動(dòng)的趨勢(shì)；同時(shí)當(dāng)慣量松弛取0.03 和0.05 時(shí)，流線圖與MRT-LBM 模型的流線圖較為接近.

為了量化以上結(jié)果，筆者測(cè)試了方腔中央的一級(jí)渦以及左右下角附近的兩個(gè)渦的渦心坐標(biāo)，結(jié)果列于表3，同時(shí)表3 也列出了不同雷諾數(shù)下有無(wú)慣量松弛因子的計(jì)算結(jié)果與基準(zhǔn)解的比較，其中a、b、c 分別代表文獻(xiàn)[11–13]中的模擬結(jié)果，從表中可以看出：不加慣量松弛因子的模擬結(jié)果與其他文獻(xiàn)的模擬結(jié)果吻合得很好，從而也進(jìn)一步驗(yàn)證了該算法程序的準(zhǔn)確性；雖然隨著慣量松弛因子R的增大，三個(gè)渦心坐標(biāo)偏移也在增大，但誤差基本都在可接受范圍以內(nèi)，當(dāng)慣量松弛因子R取0.03 到0.05 之間時(shí)，一級(jí)渦渦心坐標(biāo)誤差在5%以內(nèi)；當(dāng)慣量松弛因子R取0.1 時(shí)，一級(jí)渦渦心坐標(biāo)誤差在8%之內(nèi)；當(dāng)慣量松弛因子R取0.15 時(shí)，一級(jí)渦渦心坐標(biāo)誤差在12%以內(nèi).

圖7 不同慣量松弛因子下頂蓋驅(qū)動(dòng)流的流函數(shù)等值線圖

表3 頂蓋驅(qū)動(dòng)流的渦的位置

續(xù)表3 頂蓋驅(qū)動(dòng)流的渦的位置

4.3 慣量松弛因子在D3Q15 的LBGK 模型中的應(yīng)用

筆者還采用了LBGK 的D3Q15 模型對(duì)三維頂蓋驅(qū)動(dòng)流進(jìn)行了模擬，并將慣量松弛因子與格子Boltzmann 方法結(jié)合，研究慣量松弛因子在三維格子Boltzmann 方法中是否也具有通用性.對(duì)于三維LBGK 模型，當(dāng)雷諾數(shù)Re 取1000、ε取1.0E?12、慣量松弛因子分別取0、0.03 和0.1 時(shí)的收斂步數(shù)，如表4 所示，隨著R值增大，提高的計(jì)算效率也越來(lái)越高.

表4 不同慣量松弛因子下的收斂步數(shù)

圖8 至圖10 為不同慣量松弛因子下頂蓋驅(qū)動(dòng)流流線的不同視角圖，將不加慣量松弛因子的流線分布與文獻(xiàn)[14]中的流線分布進(jìn)行了對(duì)比，可以看出兩者不同視角的流線圖基本一致，驗(yàn)證了程序的正確性；當(dāng)慣量松弛因子取0.03 和0.1 時(shí)，空腔流線與不加慣量松弛因子時(shí)的空腔流線基本一致；但當(dāng)R值增大，其在計(jì)算效率提升的同時(shí)，中心渦的渦心坐標(biāo)也有向右上方偏移的趨勢(shì).綜上慣量松弛因子在三維格子Boltzmann 方法中一樣具有通用性.

5 結(jié)論

慣量松弛因子一般在N-S 方程中發(fā)揮作用，但在格子Boltzmann 方法中依然可以發(fā)揮極大的作用.以二維頂蓋驅(qū)動(dòng)流為例，分別采用不同雷諾數(shù)、不同慣量松弛因子進(jìn)行了數(shù)值模擬，并將計(jì)算結(jié)果與基準(zhǔn)解進(jìn)行了比較，可以看出慣量松馳因子R存在一個(gè)合理的最佳區(qū)間，當(dāng)慣量松弛因子取0.03 時(shí)，計(jì)算誤差在10%以內(nèi)；慣量松弛因子取0.05 時(shí)，計(jì)算誤差在15%以內(nèi)；當(dāng)慣量松弛因子取0.1 時(shí)，計(jì)算誤差在25%以內(nèi)；且隨著慣量松弛因子的增大，計(jì)算效率也在大幅度提高；但研究也發(fā)現(xiàn)當(dāng)慣量松弛因子取0.01 到0.02 之間時(shí)，程序收斂速度會(huì)不穩(wěn)定.綜上所述慣量松弛因子R的下限建議取0.02，并且不宜大于0.1，同時(shí)其最佳取值區(qū)間為0.03 到0.05.

筆者首次討論了將慣量松弛因子與格子Boltzmann 方法結(jié)合，研究表明它能有效解決格子Boltzmann 方法中LBGK 模型效率低的問(wèn)題，并且具有通用性，MRT-LBM 模型較LBGK 模型而言數(shù)值更為穩(wěn)定，在LBGK 模型中加入慣量松弛因子與MRT-LBM 模型所得結(jié)果較為接近，因此將慣量松弛因子引入到格子Boltzmann 方法中能在保證一定精度的前提下提高收斂速度、加快程序計(jì)算效率等，使其能在工程材料與能源環(huán)境領(lǐng)域發(fā)揮重要應(yīng)用.

圖8 R=0 時(shí)流線分布的不同視角

圖9 R=0.03 時(shí)流線分布的不同視角

圖10 R=0.1 時(shí)流線分布的不同視角