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

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

1 引言

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

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

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

2.1 慣量松弛因子與LBGK 模型

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

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

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

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

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

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

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

在上式中引入慣量項Rfk，在迭代過程中，上式的迭代格式成為

其中R為慣量松弛因子.

2) 遷移步驟如下

2.2 多松弛格子Boltzmann 模型

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

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

其中

松弛矩陣

3 物理模型

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

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

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

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

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

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

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

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

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

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

圖3 兩種模型中軸線上的速度對比

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

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

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

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

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

圖5 水平中軸線上的速度對比

圖6 垂直中軸線上的速度對比

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

圖8 R=0 時流線分布的不同視角

圖9 R=0.03 時流線分布的不同視角

圖10 R=0.1 時流線分布的不同視角