慣量松弛因子在格子Boltzmann 方法中的應用

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

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

1 引言

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

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

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

2.1 慣量松弛因子與LBGK 模型

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

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

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

平衡態分布函數為

通過Chapman-Enskog 展開，可以得到運動粘性系數與松弛頻率的關系為

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

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

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

其中R為慣量松弛因子.

2) 遷移步驟如下

2.2 多松弛格子Boltzmann 模型

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

其中f 是一個矢量，代表格子節點上的速度分布函數，m 代表一系列相互獨立的矩，meq是矩m 的平衡態值，M是一個正交的轉換矩陣，可把f 轉換為m

其中

松弛矩陣

3 物理模型

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

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

圖1 二維頂蓋驅動方腔流示意圖

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

4 計算結果及分析

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

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

表1 數值模擬結果對比

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

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

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

4.2 慣量松弛因子在D2Q9 的LBGK 模型中的應用

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

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

圖4 收斂步數與慣量松弛因子的關系曲線

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

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

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

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

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

圖7 不同慣量松弛因子下頂蓋驅動流的流函數等值線圖

表3 頂蓋驅動流的渦的位置

續表3 頂蓋驅動流的渦的位置

4.3 慣量松弛因子在D3Q15 的LBGK 模型中的應用

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

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

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

5 結論

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

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

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

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

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