多層速度格子Boltzmann 方法進展及展望

楊 鯤，單肖文

(1. 南方科技大學 前沿與交叉科學研究院，深圳 518055；2. 南方科技大學 工學院 力學與航空航天系，深圳 518055)

0 引 言

誕生于格子氣自動機[1]的格子Boltzmann 方法（LBM），與傳統直接求解宏觀量的計算流體力學（CFD）方法相比，由于避開了Navier-Stokes 方程非線性項的處理，且沒有求解壓力泊松方程這一繁重過程，使得LBM 具有程序編制簡單、并行計算效率高的優點[2]，在湍流、多相流、多孔介質流動及滲流、顆粒流、微流動、氣動聲學等領域的模擬中得到了廣泛應用[3-6]。早期的LBM 離散模型都是基于半經驗理論推導，其約束條件往往局限于滿足質量和動量守恒方程，并以低馬赫數假設為基礎，能量方程無法正確還原，使得LBM 方法只適用于等溫低速弱可壓流動成為了較普遍的觀點[2，7]。

隨著航空航天科技發展，高速流動的數值模擬在飛行器設計驗證中起到越來越重要的作用[8]，CFD 領域涌現了多種可壓縮流動的處理格式與方法，例如黎曼問題求解器和TVD 格式[9]、用于高精度激波捕捉的WENO 格式[10]，以及計算激波-湍流相互作用的Compact-WENO 混合格式[11]；這些數值方法在理論研究和工程應用領域都取得了可觀的成果[12-13]。不可否認的是，盡管傳統CFD 方法已經過了多年的發展和完善，在具體問題中要同時準確解析流動小尺度結構和激波間斷，實現模擬方法的精度、耗散、穩定性和計算效率等因素的最佳平衡，仍然存在著諸多挑戰[14]。另一方面，工程傳熱流動問題中，高效換熱器普遍具有形狀復雜的管腔結構，處理復雜幾何邊界、并兼顧精度和計算效率也對流動模擬方法提出了更高的要求。

LBM 源自于氣體動理學理論，對氣體運動描述的物理背景清晰，且在處理復雜幾何外形、特別是多孔結構時具有獨特優勢，理論上LBM 在模擬可壓縮流動和復雜結構的傳熱流動問題方面擁有巨大潛力。研究者們對傳統LBM 模型進行了諸多改進，以期突破LBM 在計算可壓縮流動和熱力學流動時的局限性，其基本思想是增加LBM 離散模型的自由度和約束條件以還原能量守恒方程。改進方案大致分為兩類，一類是額外增加一個能量（內能或總能）分布函數以滿足能量方程約束[15-20]，稱之為雙分布函數模型。然而，根據氣體動理學理論，氣體的密度、動量和能量都是氣體分子熱運動的統計結果，其動力學狀態也應當只需由單一的分布函數描述（對于單組分、單原子氣體而言），因此引入能量的分布函數必然需要與密度分布函數進行耦合，該耦合關系帶有一定的經驗性，在增加復雜度的同時也對模型的廣泛適用性帶來了一定的限制。另一種方案是構造更高階的單一平衡態分布函數，在平衡態分布函數中引入更多的待定系數，并使用更多的離散速度點，通過增加離散點數量（未知數個數）的方案來滿足額外的能量方程約束，典型做法是在每個離散速度方向上使用2 個以上的速度大小模態，稱之為多層速度模型。與雙分布函數模型類似，多層速度模型在能量方程約束條件、分布函數的形式和離散速度點選取上也具有一定的經驗性，自Qian[21]和Alexander 等[22]提出以來，出現了形式多樣的方案且各有優缺點。

多層速度模型按平衡態分布函數展開式的構造方法主要分為三類：第一類遵循傳統LBM 模型的思路，將Maxwell-Boltzmann 分布函數泰勒展開至更高階數[23-26]。第二類由Shan 等[27]提出，他們發展了Grad[28-29]的早期工作思路，認為平衡態分布函數可由Hermite 多項式級數展開，只需要將展開級數保留足夠的截斷階數，即可分別還原到Navier-Stokes 方程的質量、動量和能量方程，而密度、動量、能量等宏觀統計量均可由Hermite 展開式的各階系數組合表達，并在進一步的工作中闡明傳統等溫弱可壓LBM是Hermite 多項式模型的低階形態，且離散速度是五階精度Gauss-Hermite 積分公式的積分點的事實[30]。第三類由Xu 等[31]提出，其基本思想是不預設平衡態分布函數的展開形式，而是將離散模型還原質量、動量和能量方程的約束條件視為平衡態分布函數離散點的線性代數方程組，直接求解方程組得出離散分布函數。

由于Maxwell-Boltzmann 分布函數的各階Hermite多項式展開是確定的，Gauss-Hermite 積分公式的積分點也具有確定性，使得上述第二類基于Hermite 多項式構造的多層速度模型不依賴于經驗而具有普適性。然而，由于該模型的數學基礎理論較為晦澀，各部分細節的闡述散落于跨度接近20 年的多篇文獻中，使得科研工作者對該模型進行深入理解存在一定的困難。為了彌補這個缺點，使讀者連貫、系統地理解Hermite 多項式模型的構造原理和方法，同時也在與其他LBM 模型的對比中全面認識各類方案的優劣，我們對Hermite 多項式模型重新梳理歸納的同時，也一并列出了其他多層速度模型的構造思想，作為LBM 多層速度模型的一個簡要總結性工作。

本文的各部分內容為：第1 節為Boltzmann 方程和LBM 的基礎理論。由于多層速度模型起源于經典的等溫弱可壓LBM，且經典LBM 模型的很多思想在多層速度模型中仍然適用，故第2 節對等溫弱可壓LBM 進行簡單介紹。第3 節介紹各類多層速度LBM模型（除Hermite 多項式模型以外）的基本思想。第4 節詳細描述Hermite 多項式模型的理論和平衡態分布函數的具體形式，以及基于Gauss-Hermite 積分公式的離散速度模型及其構造方法。第5 節對LBM 和傳統CFD 方法的結合進行了簡要介紹，例如LBM 有限差分、LBM 有限體積和LBM 有限元方法。第6 節對現有的多層速度模型存在的問題進行總結，并提出未來需要完善和發展的方向。

1 Boltzmann-BGK 方程

由于無量綱化的Boltzmann-BGK 方程和平衡態分布函數分別與式（10）和式（11）形式上完全一致。因此

2 等溫弱可壓LBM 模型

等溫弱可壓LBM 模型是最早提出的LBM 模型，并在不可壓流動計算中得到廣泛應用。嚴格意義上，包括氣體、液體在內的所有流體都具有可壓縮性，例如液體的壓縮性可用Tait 方程描述[34]。傳統求解壓力泊松方程的不可壓流計算方法將導致無限大聲速，也就是壓力擾動的無限速度傳播，這與物理事實是相悖的。采用弱可壓方法計算“不可壓流動”理論上更加接近物理本質。本節對等溫弱可壓LBM 模型進行簡要陳述，這也將成為多層速度LBM 模型修正和拓展的基礎。

2.1 速度空間離散模型

圖1 等溫弱可壓流動LBM 的一維和二維離散速度模型Fig. 1 One- and two-dimensional discrete velocity models for isothermal weakly compressible flows

2.2 時間和空間離散

3 多層速度LBM 模型

從第2 節可以看到，等溫弱可壓流動LBM 模型，其離散速度在每個方向都只有一個值（分布函數信息只傳播到相鄰的網格點），因此稱之為單層速度模型。推導該模型時，只約束了離散模型的質量、動量和動量輸運量，故無法準確還原能量方程。受此啟發，業內學者開始對離散模型施加能量相關（f(0)的ξ三階以上矩）的約束條件，伴隨而來的是更高階的f(0)展開式和更多的離散速度點，其中普遍做法是在各離散速度方向上采用2 個以上的速度值（模態），因此這些模型稱之為多層速度模型。本節對幾類較為典型的多層速度模型進行簡要陳述。

3.1 早期多層速度模型

圖2 Qian 的二維離散速度模型[21]Fig. 2 Two-dimensional discrete velocity model of Qian[21]

圖3 Alexander 等[22]的二維離散速度模型Fig. 3 Two-dimensional discrete velocity model of Alexander et al[22]

3.2 Watari-Tsutahara 模型

3.3 比熱可變模型

圖4 Chen 等的二維多層速度模型[23]Fig. 4 Two-dimensional discrete velocity model of Chen et al.[23]

圖5 Watari-Tsutahara 二維多層速度模型[24-25]Fig. 5 Two-dimensional discrete velocity model of Watari-Tsutahara[24-25]

圖6 Kataoka-Tsutahara[26]二維多層速度模型，圖中的坐標代表離散速度點ξ a/cs在第一象限的坐標Fig. 6 Two-dimensional discrete velocity model of Kataoka-Tsutahara[26]. Coordinates represent the discrete velocities in the first quadrant

圖7 離散Boltzmann 方法的幾種二維離散速度模型[39]Fig. 7 Several two-dimensional discrete velocity models in DBM[39]

4 Hermite 多項式模型

4.1 速度空間展開

從第3 節的介紹可以看出，大部分多層速度模型還是沿襲了單層速度模型“預設平衡態分布函數離散形式—預設離散速度模型—待定系數法求解”的思路。這種思路可以針對每一類特定問題都構造相對優化的多層速度模型，但經驗性較強，對于開始接觸多層速度模型的學者而言，面對種類繁多的離散模型、平衡態分布函數形式和約束條件，如何選擇最優方案往往無所適從；即便獲得了對某一問題的優化組合，當流動參數改變后，計算結果不一定同樣令人滿意，無法保證模型的普適性。

Shan 等[27,30]在Grad[28-29]工作的基礎上，發展了Hermite 多項式模型，該模型基于Hermite 多項式級數展開的數學性質，使平衡態分布函數的展開形式和離散速度模型都不再依賴于經驗而具有確定性，且可以證明單層速度模型事實上是多層速度模型的低階形態。這正是本節將詳述的內容。

參照Grad 的定義，D維空間的n階Hermite 多項式定義為n階張量：

表1 不同流動分布函數的Hermite 截斷式最低階數Table 1 Lowest order of truncated Hermite distribution function for different types of flows

4.2 時間和空間離散

4.3 離散速度模型

圖8 2 種D2V17 離散速度模型Fig. 8 Two types of D2V17 discrete velocity models

在之前的討論中，有一個重要的問題沒有展開——如何獲得具有K階精度的數值積分公式（式點序列要保留表2 中列出的至少2 個8 積分點組，其總積分點數最少為 1+4×5+8×2=37。若保留表2中的前8 組積分點組與權重，則得出與Pillippi[51]相同的D2V37 模型，如圖9 所示。

表2 二維積分公式的積分點與權重Table 2 Quadrature abscissas and weights of the 2D integral formula

圖9 D2V37 離散速度模型Fig. 9 D2V37 discrete velocity models

對于三維情況下的積分公式，也可以遵循同樣的思路求解。相關的細節過程可參考文獻[47]。

4.4 等溫弱可壓模型的推導

5 與傳統CFD 方法的結合

經典的LBM 中時間和空間離散均為沿特征線積分（參考2.2 節和4.2 節），并設置歸一化的 Δt=1，這極大簡化了LBM 的計算程序。然而這種時空離散方法也帶來了幾個弊端：

1）由于 Δx為固定值，只能采用均勻正方形（體）網格。對于邊界層流動計算而言，為了解析邊界層內的流動而加密全計算域網格會極大增加計算量。在計算量和解析度之間取平衡的辦法是使用局部加密網格，但是粗網格和細網格之間的插值會帶來新誤差[52]。

2）無法使用貼體網格。在傾斜方向邊界，或者曲面邊界上，需要同時使用笛卡爾網格和浸沒邊界法[53]，以鋸齒形邊界近似斜邊界或者曲面邊界。這種非貼體網格由于精度較低，且在對流擴散問題中易造成收斂性困難[54]和穩定性問題，往往需要與局部加密網格和多松弛時間模型（MRT）結合[55]。

3）計算穩定性。經驗表明，經典LBM 方法在雷諾數較大時只能使用較小的物理時間步長才能獲得穩定的數值解；當使用Zhou-He 非平衡態反彈邊界條件[56]時穩定性問題更加突出，這對實際工程湍流的計算會形成較大的障礙。

由于傳統CFD 方法可以使用邊界層加密貼體網格，使用合適的時間與空間格式時，數值穩定性良好，將傳統CFD 方法與LBM 結合成為了解決以上問題的一種思路。Reider 和Sterling[57]、Cao 等[58]將有限差分法與LBM 模型結合，提出了FDLBM 方法。FDLBM 的基本思路是將LBM 方程（式（104））視為關于各個離散分布函數fa的雙曲型方程，通過有限差分法的時間離散格式，例如Euler 格式或者Runger-Kutta 格式，將 ?fa/?t離散化；并使用迎風型或中心型格式將對流項 ξa·?fa進行空間離散[59]。在可壓縮問題中，需要將多層速度模型和激波捕捉差分格式結合。Kataoka 和Tsutahara[26]在多層速度模型的基礎上，使用二階迎風格式離散對流項；Wang 等[60]將二階TVD 和五階WENO 格式[10]與Kataoka-Tsutahara多層速度模型結合使用；許愛國等[39]則使用了NND 格式[61]來處理對流項，并在一維和二維黎曼問題的計算中取得較為滿意的結果。

使用隱式時間離散是傳統CFD 方法中提高數值穩定性的常用方法，伴隨而來的是需要迭代求解線型代數方程組，這往往使得計算量和實施難度顯著增加。在FDLBM 中，Guo 和Zhao[62]使用了與經典LBM 類似的隱式時間處理辦法，即類似于式（30）定義一個中間時刻的fˉa代入FDLBM 方程中，可將隱式方法轉換為顯式方法。Wang 等[60]沿用類似的思路，將隱式-顯式Runger-Kutta 格式[63]引入FDLBM 并消除了隱式格式需要的迭代求解過程。

為了適應復雜幾何構型和不規則形狀網格的計算需求，Nannell 和Succi[64]、Benzi 等[65]將有限體積方法引入LBM，提出了FVLBM 方法。FVLBM 的基本思想是，將LBM 方程（式（104））對流場中任意一個控制體單元 Ω （邊界為 ?Ω ，此處 Ω不代表碰撞模型）進行體積分得到：

則式（123）可寫為關于分布函數單元體平均量和單元面對流通量的方程。求解該方程的過程即與傳統有限體積方法類似。Chen[66]和Peng 等[67]將FVLBM的理論和其在非均勻網格上的應用進行了完善。為了提高穩定性，Stiebler 等[68]使用迎風型格式重構單元面通量。Patil 和Lakshmisha[69]則將TVD 格式應用于單元面通量的重構。

隨著計算流體力學有限元方法日趨成熟[70-71]，Lee 和Lin[72]、Shi 等[73]將有限元方法引入LBM，發展了 有 限 元LBM 方 法。近 年 來，MacMeccan 等[74]、Krüger 等[75]使用混合型有限元LBM 方法模擬了可變形顆粒流動并應用于生物流計算領域。由于主題和篇幅所限，此處不對有限元LBM 方法作進一步展開。有興趣的讀者可參閱有關文獻了解相關細節。

6 總結與展望

LBM 多層速度模型起源于單層速度模型，在熱流和可壓縮流動應用需求下不斷發展，經歷了平衡態分布函數展開從二階提高到四階、離散速度模型從單層網格提高到2 層或者更多層網格、約束條件也逐步增多的過程。Hermite 多項式模型則提供了一個既理論又實用的判定框架：要想準確還原能量方程，平衡態分布函數必須展開到速度四階量。Watari-Tsutahara預設的多層速度模型[24-25]與該結論殊途同歸。Hermite 多項式模型也同時指出，離散速度模型事實上為相應階數積分公式的積分點和權重，只需查表（或用程序包計算）找出這些積分點和權重，選出一組應用即可，而不是根據經驗預設和優化。盡管Hermite 多項式模型給出了較為簡潔和確定的形式，但是依照經驗理論進行馬赫展開或待定系數得到的其他多層速度模型，在特定范圍和領域內仍具有一定的應用價值。

必須承認，相對于計算等溫弱可壓流動的經典單層速度模型而言，多層速度模型仍然是一種尚不成熟的模型，在很多方面仍然需要改進和完善。這主要在以下幾個方面：

1）模型的誤差、色散與耗散、穩定性分析。目前對單層速度模型的誤差和穩定性分析已有不少工作，例如：Sterling 和Chen[76]對D2Q7（六邊形模型）、D2Q9 和D3Q15 進行線性穩定性分析后得出 τ ＞0.5才能穩定的結論。Marie 等[77]對D3Q19 模型進行色散和耗散分析，并與低色散高階差分格式進行了比較。Silva 等[78]和Bauer 等[79]對最常用的D3Q15、D3Q19和D3Q27 模型進行了截斷誤差分析，指出這三種模型都能在低馬赫數下還原動量方程，但是動量輸運項的高階截斷誤差不同，這將導致在各向異性流動的計算中出現較大差別；在將平衡態分布函數展開式根據連續平衡態分布函數的二階矩約束條件和D3Q19 速度模型進行系數修正后，能獲得與D3Q27 模型相當的各向同性水平。

對多層速度模型而言，McNamara 等[80]對D3Q21模 型（D3Q15 加 上 (±2c,0,0)、 (0,±2c,0) 和 (0,0,±2c)這6 個離散點）進行了時間離散格式分析，指出使用Lax-Wendroff 時間格式將在高瑞利數下具有更好的穩定性。Wissocq 等[81]則對D2Q9 和D2V17 進行了模態和穩定性分析，指出D2V17 能準確解析線性等溫Navier-Stokes 方程的3 個特征值模態，而D2Q9 則在這些模態的高波數上產生一定的誤差。這一點與Hermite 多項式模型的結論是一致的，即只有使用七階積分公式的積分點（即D2V17 模型）作為離散速度點才能準確還原動量方程；否則由于速度三階量誤差的存在而必須使用低馬赫數假設。

值得注意的是，在Hermite 多項式模型框架下，經典D3Q15、D3Q19 和D3Q27 模型從積分精度角度來講，并不存在差別，因為都具有五階積分精度；從計算經濟性原則來講，應該選擇積分點最少的D3Q15 模型。然而Silva 等[78]和Bauer 等[79]的工作提示我們，在選擇離散模型時還需要考慮各向同性、截斷誤差等更多方面的因素，即計算最經濟的模型不一定是最優的計算結果。Sieber 等[82]將Hermite 多項式模型與D2Q9 模型、Chen 等[23]的多層速度模型進行了線性穩定性分析對比，發現對等溫弱可壓流動而言，將平衡態分布函數展開到三階Hermite 多項式級數，或使用D2V17、D2V37 模型都能顯著提高穩定性；對于可壓縮熱流而言，使用四階Hermite 展開的f(0)與D2V37 模型，穩定性都優于Chen 等[23]提出的多層速度模型。該結論在Hermite 多項式模型的理論框架下是易于被推斷的。然而，對于完全還原了能量方程、且離散速度點數量最少的D2V37、D3V103 模型[47]而言，與同階積分精度、但具有更多離散速度點的模型橫向對比，例如D3V107 模型，是否也存在各項同性、穩定性、色散和耗散差異？這一方面需要在后續工作中展開詳盡地理論分析；另一方面也需要大量的數值試驗，例如方管Poiseuille 流動、旋轉管道流、強可壓縮流動、大溫差熱流等實際流動計算，來驗證這些多層速度模型之間的性質差異。

2）邊界條件。LBM 單層速度模型在經過多年發展后，其邊界條件處理方法已經日趨完善，有關工作可謂汗牛充棟[33,36,83]，對于壁面反彈邊界條件[56,84]、周期邊界條件[85]、適用于流動出入口的無反射邊界條件[86-87]以及其他類型的邊界條件等都已有充分的研究工作，可滿足絕大多數計算工況的需求。然而，當使用多層速度模型時，一方面由于離散速度點跨越多層網格，使得邊界附近的內部結點與邊界結點區分變得模糊；另一方面平衡態分布函數的形式也變得更為復雜，如何恰當處理邊界條件成為了一個難點，相關的工作也屈指可數[83]。Malaspinas 等提出[88]，將邊界節點的分布函數離散點按內部和邊界分為已知量和未知量，將離散分布函數各階矩的約束條件構成方程組，通過求解方程組構造一種“通用”的邊界條件處理方法。他們在D2Q9 和D3Q19 模型中進行了測試，并指出該方法可適用于多層速度模型的邊界條件（但沒有后續驗證）。Meng 和Zhang[89]提出了適用于多層速度模型的擴散-反彈邊界條件格式，并應用到D2V16、D2V17、D2V37 和D3V121[45]這四種模型，對等溫流動和熱流進行了測試計算。Lee 等[90]將Hecht 和Harting 針對D3Q19 的非平衡態壁面反彈邊界條件[91]，推廣到了用多層速度模型計算等溫弱可壓流動的情況，并用D2V17 模型進行了驗算。Klass等[92]進一步延續了Lee 等的工作，并將其應用到了弱可壓熱流中，使用D2V17 和D2V37 模型進行了驗證，計算結果比Meng 和Zhang 的擴散-反彈邊界條件表現更優。總體來看，針對多層速度模型的邊界條件研究工作還處于起步階段，且大部分局限于弱可壓熱流和Dirichlet 邊界類型；對于壓縮性較強的跨聲速、超聲速流動，以及流動出入口涉及的無反射邊界條件，尚無相關工作涉足。將單層速度模型中的各類邊界條件處理方法擴展到多層速度模型，并在各類流動中廣泛驗證測試，仍是今后需要較多投入的方向。

3）高性能算法。LBM 最大優勢是計算速度和高度并行特性，然而它也有一個眾所周知的缺點：由于需要存儲數量眾多的離散分布函數，對計算機內存的使用量比傳統CFD 方法更多，因此有不少工作聚焦于如何改進算法，減少內存消耗并提高計算效率[93-95]。另外，正如上一節所述，經典的LBM（相對于FDLBM 和FVLBM 及其他混合方法而言）一般使用均勻網格，當涉及到邊界層計算時，將LBM 與自適應局部加密網格（AMR）結合[96-98]是更加經濟有效的方案。這些算法與目前流行的GPU 高性能計算相結合，能發揮LBM 的最大計算潛力[99-101]。然而，目前這些高性能算法都是集中應用在單層速度模型中。對于多層速度模型而言，特別是三維情況，準確模擬可壓縮熱流的最經濟模型是D3V103，但其分布函數離散點的數量已經是D3Q27 的4 倍。這一方面使其內存占用問題更加突出；另一方面，在并行化和局部自適應加密網格時，不同塊之間的界面數據交互量也更大，需要相當謹慎地處理。然而這些往往是工程實際應用中決定算法效率的制約因素。發展適用于多層速度模型的高性能算法，包括內存節約算法、高效率通信和AMR 方法、CPU+GPU 異構高效算法，并在各類復雜幾何外形、復雜流動中廣泛測試，應引起業內學者的廣泛關注。

總體來說，LBM 模型，特別是多層速度模型，仍然需要諸多補充和完善才能成為一種可靠、高效的流體力學研究手段和應用技術。正如1998 年Chen 和Doolen[2]在文章結尾中提到，LBM 多相流和復雜反應系統模型主要集中于等溫弱可壓流問題，適應于可壓縮熱流的LBM 模型仍待開發。20 余年后的今天，這些方面雖然取得了一些進步，但是離可靠性要求還有不少差距。我們期望多層速度模型能夠百花齊放，特別是具有較完備數學理論框架的Hermite 多項式模型，能在這些方面繼續發展壯大并做出應有的貢獻。

附錄A