高速可壓反應流的二維簡化離散玻爾茲曼模型

林傳棟

(中山大學 中法核工程與技術學院, 珠海 519082)

0 引 言

化學反應流不僅廣泛存在于自然界中，更在人類日常生活、工業生產、航空航天和國防建設等諸多領域發揮著越來越重要的作用[1-3]。在化學反應流中，各類流體動力學行為和化學反應過程相互耦合，具有典型的流體力學非平衡、熱力學非平衡和化學非平衡效應，表現出豐富多變的物理化學現象[1-3]。對于此類非定常、非平衡的復雜物理化學系統，不能使用定常流場的解析方法或數值模型進行研究，而需要使用包含多物理場耦合和時空演化過程的物理模型。尤其是高超聲速飛行器[2-4]還同時面臨稀薄氣體效應、化學反應效應和壁面熱流效應等外界因素的影響。這些極端復雜條件對如何有效模擬和預測高超聲速飛行器的運行過程提出挑戰。因此，高速反應流的數值建模與模擬研究在航空航天等領域具有特別重要的實際應用背景，是進一步提升航空航天工程技術的關鍵問題之一。

相對于高超聲速飛行器，航天飛行器再入過程跨越了更加寬泛的時空尺度。具體來講，航天飛行器以（超）高馬赫數依次經歷稀薄自由分子流區、過渡流區、滑移流區以及連續流區，這一經歷伴隨著復雜的變密度、變溫度、變速度、電離和化學反應等真實飛行環境。在各個流體區域，氣體物性表現出顯著差異，其繞流狀態具有典型的跨尺度效應。同時，幾何尺寸差異明顯的外形部件也會使飛行器不同區域表現出多流區共存的多尺度混合流動現象。在高超聲速再入過程中，還具有強黏性、高摩阻、高熵層、強氣動加熱、熱力學與化學非平衡相互耦合等典型特征[5-6]。為此，需要構建準確高效、穩定可靠的跨尺度模型，用以研究分析高超聲速飛行器跨流域飛行的氣動力熱機理，明晰復雜繞流的時空變化特點以及非平衡熱化學行為表征。

然而，傳統的宏觀數值研究大都基于連續性假設，側重于模擬研究系統宏觀演化過程，忽略了實時伴隨流體宏觀行為的豐富的熱力學非平衡效應，導致其應用范圍不夠廣泛、數值仿真結果不夠精確。目前，基于Boltzmann 方程的動理學模型具有微介觀尺度的優勢，相關方法的應用研究已經成為國際國內的熱點課題，包括格子氣自動機（lattice gas automaton，LGA）、格子玻爾茲曼方法（lattice Boltzmann method，LBM） 、 離 散 玻 爾 茲 曼 方 法（ discrete Boltzmann method，DBM）、離散速度方法、離散統一氣體動理論方法、氣體動理論統一算法和矩方法等。

作為一種元胞自動機，LGA 不僅是解方程（恢復方程）的一種全新的數值方法，更是研究各類復雜系統的一種基本工具[7]。最初，LGA 由法國的Hardy、Pomeau、De Pazzis 于1973 年首次提出，目標是能夠模擬Navier-Stokes（N-S）方程組描述的流體系統，但是由于該模型旋轉不變性的缺失，導致它具有高度各向異性[8]。后來，Frisch、Hasslacher、Pomeau 在1986 年優化了LGA 模型，其格子矢量具有四階張量對稱性[9]。時至今日，LGA 已經被成功應用于流體力學、化學、電磁學、熱聲學等領域。

LGA 的基本思想是：不同的微觀行為可以導致相似的宏觀表征，虛擬粒子在網格上遷移和碰撞的過程保持物理量守恒，如質量、動量、能量守恒[7]。為了提高計算效率，這種虛擬世界的微觀動力學應該盡可能簡單[7]。此外，LGA 可看作是LBM 和DBM 的鼻祖。廣義上講，LGA 及其后代動力學模型被認為具有許多相似優點：1）設計方案簡單易行，因此便于編寫代碼和檢驗程序。2）碰撞規則具有局部性，因此非常適合大規模并行計算。3）可通過選取合適的邊界條件來模擬具有復雜幾何結構的系統。然而，LGA具有統計噪聲等問題，而LBM 和DBM 能夠有效避免該問題。

近30 年來，LBM 被成功應用到諸多復雜流體系統[10-12]，并且在早期就被拓展到化學反應流領域。早在1997 年，歐洲院士Succi 等[13]首次提出燃燒系統的LBM，其限制條件包括：化學反應極快、流場不可壓、反應熱對流場沒有影響。2015 年，日本的Yamamoto 等[14]使用LBM 模擬了三維空隙結構系統中的燃燒現象。2018 年，德國Hosseini 等的課題組[15]構建了多組分系統的LBM，并解決了擴散方程誤差引起的質量不守恒問題。同時，國內學者也在相關領域做出了突出貢獻。華中科技大學的陳勝等[16]把LBM 應用到二維和三維低馬赫數燃燒系統，將流體密度和溫度場直接耦合。西安交通大學的陶文銓院士、何雅玲院士等的課題組[17-18]使用LBM 研究包含多相、多組分、化學反應的能源轉換系統。另外，國內外還有許多其他優秀科研團隊將LBM 應用到包含溶解、沉淀、流固耦合等復雜因素的反應流系統[19-21]。

盡管傳統LBM 的應用研究已經取得了巨大進步，但基本都是還原傳統宏觀控制方程的解法器，基本沒有聚焦于詳細的熱力學非平衡信息，也沒有實現真正意義上的化學反應和可壓縮流體的耦合。作為一種新型的微介觀動理學方法，近幾年發展起來的DBM 能夠解決物理精度與計算效率之間的矛盾。DBM 是玻爾茲曼方程在速度空間的特殊離散化形式，可以看作傳統LBM 的進一步發展或變異[22-36]。值得說明的是，DBM 繼承了玻爾茲曼方程描述非平衡物理系統的主要功能。可以證明，一方面DBM 在連續性極限條件下能夠恢復宏觀流體控制方程組；另一方面DBM 還包含宏觀流動方程組之外的熱力學非平衡效應[22-36]。因此，DBM 可用于處理宏觀模型失效、非平衡效應較為顯著的微尺度物理系統、稀薄氣體系統和滑移層流體區域等相關問題[31]。

目前，DBM 在化學反應流方面已經得到初步應用。2013 年，DBM 已被初步發展到爆轟領域，并被用于模擬分析爆轟波附近的非平衡效應[22]。2014 年，林傳棟等[23]構建并使用極坐標系DBM 模擬研究了內爆和外爆過程的非平衡強度，并從動理學角度分析了粒子速度分布函數的主要特征。2015 年，許愛國等[24]構建并使用多松弛DBM 模擬研究了黏性和熱傳導在爆轟系統中對局域和全域非平衡效應的影響。2016 年，張玉東等[26]使用DBM 模擬研究了負溫度系數對爆轟系統的影響機理。之后，林傳棟等構建了多組分DBM 用于模擬預混、非預混和部分預混反應流系統[27]，又使用高效、準確的多松弛DBM研究了不穩定反應熱流演化過程的非平衡行為[30]。

為了在連續性極限條件下恢復N-S 方程組或者更高階的流體力學方程組，并且為了包含更準確和詳實的熱力學非平衡信息，要求DBM 包含足夠多的動理學矩關系。現有的DBM 模型至少能夠恢復N-S方程組，與之對應，至少包含了7 個張量形式的動理學矩關系。例如，對于二維DBM，當包含額外自由度時，需要構建16 速度的離散速度模型[30]，與之對應，需要包含16 組形式統一的DBM 演化方程。然而，在實際的物理系統中，有時（或有些區域）物理場蘊含的非平衡效應可能較弱，不需要特精確和太費時的模擬工具，而需要使用較為簡潔和高效的模型。

為此，本文首次提出并構建了簡化的二維DBM。該模型只包含9 個獨立的動理學矩關系，并且在連續性極限條件下能夠恢復帶有化學反應的Euler 方程組。該DBM 包含9 個離散速度，基于9 組形式統一的DBM 演化方程。相對于之前構建的16 速度模型[30]或者24 速度模型[24]，該DBM 具有更高的運算效率。本項目的開展有助于推動動理學方法在化學反應流的發展應用，并進一步促進其在復雜多尺度物理系統的模擬研究。

1 二維簡化DBM

1.1 離散玻爾茲曼方程

圖1 離散速度示意圖Fig. 1 Sketches for the discrete velocity model

1.2 動理學矩關系

1.3 化學反應進程

為了簡單起見，本文忽略電離、熱輻射、活性自由基和可逆反應等效應[25]；同時，考慮較為一般的情況，即化學反應的時間尺度大于熱力學弛豫時間，且小于流體系統的特征尺度[25]。那么，由此可以進一步假設：化學反應引起的粒子速度分布函數的改變率近似等于平衡態速度分布函數的改變率，且在化學反應的時間尺度內，局部流體的密度和速度尚未來得及改變[25]。

1.4 流體力學模型和非平衡物理量

可以看出，在物理描述功能方面，簡化二維DBM除了具有Euler 模型描述流體系統演化的功能之外，還具有描述一定熱力學非平衡行為的功能。因此，該DBM 動理學模型具有超越Euler 模型物理描述范圍的優勢。也就是說，在該簡化二維DBM 中，DBM等價于帶有化學反應的Euler 方程組和包含少量非平衡效應的粗粒化模型。

需要說明的是，對于二維DBM，當對應包含額外自由度的N-S 方程組時，需要構建16 速度[30]（或24 速度[24]）的離散速度模型，與之對應地需要包含16 組[30]（或24 組[24]）形式統一的DBM 演化方程；當對應Burnett 方程組時，至少需要26 個離散速度和26 組演化方程[36]，計算量相應增大。然而，在實際的物理系統中，有時（或有些區域）物理場蘊含的非平衡效應可能較弱，不需要特精確和太費時的模擬手段，反而使用較為簡潔和高效的模型更為可行。此時，僅有9 個離散速度和9 組演化方程的簡化DBM 更加可取。總之，簡化DBM 的優勢為：離散速度數量較少、計算效率高。

另外，在已有的二維DBM 中，平衡態分布函數和反應項都至少各自滿足16 個[30]（或24 個[24]）獨立的矩關系，DBM 等價于帶有化學反應的N-S 方程組和包含較多非平衡效應的粗粒化模型。尤其對于文獻[36]，平衡態分布函數滿足26 個獨立的矩關系，該DBM 等價于Burnett 方程組和包含更多非平衡效應的粗粒化模型。在本文構建的簡化模型中，平衡態分布函數和反應項都僅滿足9 個獨立的矩關系。相比于已有的DBM，該DBM 的局限性為：滿足的動理學矩關系較少、物理精度較低，僅適用于熱力學非平衡程度較弱的物理系統。

1.5 數值離散格式

DBM 使用離散玻爾茲曼方程描述化學反應流的演化過程，該方程形式簡單、易于編程。同時，該方程具有時空局域性特點，易于實現并行化，方便在高性能計算集群運行。另外，對于離散玻爾茲曼方程式（1）的時間和空間偏微分，可以借用已有的數值格式處理。本文使用的是二階精度的Runge-Kutta 格式和二階精度的無波動、無自由參數的耗散有限差分格式（nonoscillatory and nonfree-parameter dissipation finite difference scheme，NND）[38]分別用于處理方程式（1）的時間和空間偏微分。

1.5.1 Runge-Kutta 格式

本文采用Runge-Kutta 格式處理時間離散化問題，具體計算流程見圖2。

圖2 二階Runge-Kutta 格式流程圖Fig. 2 Flow chart of the second-order Runge-Kutta scheme

1）程序的開頭是初始化，根據物理場的初始條件，輸入參數；

2）根據給定的時間步數，即程序循環周期，判斷是否終止運行；

3）根據給定的時間間隔，判斷是否輸出數值結果；

需要說明的是，上述NND 格式具有二階空間精度，是二階中心格式和二階迎風格式的結合，是具有負系數的四階耗散項。因此，該NND 格式可以抑制奇偶失聯震蕩，能夠有效捕捉沖擊波和爆轟波[38]。

2 數值驗證

本部分通過3 個典型算例對DBM 進行數值驗證。首先，通過均勻化學反應系統驗證本模型能夠將化學反應與物理場自然耦合，模型既適用于吸熱反應，也適用于放熱反應，并且反應流系統的比熱比可調；其次，通過Sod 激波管驗證本模型適用于可壓縮流體系統，并且能夠同時捕捉稀疏波、物質界面和沖擊波的時空演化；最后，通過爆轟波算例開展了網格無關性驗證，并驗證本模型適用于超聲速可壓縮化學反應流。

2.1 均勻化學反應系統

為了驗證本模型能夠將化學反應與物理場自然耦合，本節首先考慮均勻化學反應系統：在初始時刻，化學反應物均勻、靜止分布在計算區域，初始溫度為T0=2。經過一段足夠長的時間后，化學反應結束，反應物完全轉變成產物，化學能完全釋放并轉變成系統內能。由于該系統為均勻系統，每處的物理量分布都相同，所以為了提高計算速度，可以只用一個網格點，即NX×NY=1×1。 松 弛 時 間 τ=4×10?6，空間步長Δx=Δy=1×10?4，時間步長 Δt=2×10?6，離散參數va=?3.7、vb=?2、vc=?1.5 ， ηa=4、 ηb=0、 ηc=0。另外，對計算區域四周都采用周期邊界條件。

圖3 展示了化學反應后不同化學能對應的系統溫度，其關系見式（37），此處設置比熱比 γ=1.4。符號代表DBM 模擬結果，實線代表理論解：

可以看出，DBM 的模擬結果與理論解完全吻合。而且，圖3 中化學能的數值涵蓋負值到正值的范圍。其中負值區域代表化學反應吸熱，正值區域代表化學反應放熱。也就是說，該DBM 既適用于化學吸熱系統，也適用于化學放熱系統。

圖3 化學反應后不同化學能對應的系統溫度Fig. 3 Temperature after chemical reaction under various chemical energies

下面考慮固定化學能不變Q=1，調節比熱比數值的情況。圖4 展示了化學反應后不同比熱比對應的系統溫度。符號代表DBM 模擬結果，實線代表理論解（式（37））。同樣，DBM 給出的模擬結果與理論解完全吻合。因此，數值模擬證明該模型能夠適用于不同比熱比的化學反應流系統。

圖4 化學反應后不同比熱比對應的系統溫度Fig. 4 Temperature after chemical reaction under various specific heat ratios

2.2 激波管

圖5 展示了在t= 0.02 時刻激波管中的物理量分布，圖5（a～d）依次是密度、溫度、壓強和水平速度。符號代表DBM 模擬結果，實線代表理論解。可以看出：在最左側的是稀疏波（向左傳播），在中間的是物質分界面，在最右側的是沖擊波（向右傳播）。其中，稀疏波在水平方向的跨度最大，各物理量都有清晰的空間梯度；在物質界面兩側，密度和溫度具有明顯空間梯度，而壓強和速度均勻分布；在沖擊波波陣面，各物理量梯度最大，具有劇烈的空間變化。在以上演化過程中，DBM 與理論解具有很高的匹配度。這充分說明，DBM 能夠很好地捕捉稀疏波、物質界面和沖擊波的空間分布。

圖5 物理量在激波管中的分布Fig. 5 Profiles of physical quantities in the shock tube

2.3 爆轟波

爆轟是一種特殊的化學反應流現象，在其化學反應區前沿是以超聲速傳播的沖擊波[1]。在爆轟演化過程中，密度、溫度和壓強急劇上升，流場產生劇烈變化。物理量時空變化如此劇烈，在數值模擬中容易產生數值振蕩，所以對爆轟波的數值模擬方法要具有較好的魯棒性。本節將驗證我們的DBM 模型能夠有效模擬爆轟波的傳播過程。

下面考慮沿水平方向傳播的爆轟波。在一個水平放置的長度為l0=0.06的直通道內，初始物理場設置如下：

下標L表示左側區域 0 ≤x＜0.99l0，下標R表示右側區域 0.99l0≤x≤l0，左右兩側的物理量滿足Hugoniot關系。設置化學反應熱Q=2， 馬赫數Ma=2.12643，比熱比 γ=1.4 ， 松弛時間 τ=5×10?6， 離散參數va=?3.7、vb=?2、vc=?1.5， ηa=4、 ηb=0、 ηc=0。另外，上下采用周期邊界條件，左右采用出口邊界條件。

首先進行網格無關性驗證。圖6 展示了不同空間步長的密度分布圖，即 Δx1=5×10?5、 Δx2=1×10?4、Δx3=2×10?4、Δx4=4×10?4；與之對應的水平網格點數分別是NX=1 200、 600、 300、 150。可以看出，隨著空間步長的減小（網格數的增加），模擬結果之間的差異越來越小。尤其是在 Δx1和 Δx2下的模擬結果幾乎完全重合。因此，考慮到計算效率，可以使用空間步長 Δx2進行數值模擬。同樣，可以對不同的時間步長進行數值對比。結果表明，時間步長 Δt=2×10?6和4×10?6之間的模擬結果差異可以忽略。

圖6 在不同空間步長下爆轟波附近的密度分布圖Fig. 6 Profiles of density around the detonation wave under various spatial steps

圖7 展示了在空間步長Δx=1×10?4、 時間步長Δt=2×10?6下爆轟波附近的模擬結果。圖7（a～d）依次對應密度、溫度、壓強和水平速度的空間分布，其中符號代表DBM 模擬結果，實線代表Zeldovich-Neumann-Doering（ZND）理論解[1]。可以看出，DBM模擬結果與ZND 理論解具有非常好的一致性。尤其是當化學反應結束后，在爆轟波后側形成穩定的“平臺區域”，此處的DBM 模擬結果為 ρ=1.48065、T=2.06216、p=3.05334、ux=?1.69974。對比ZND 理論解，相對誤差依次為： 0 .015%、 0 .047%、0 .032%、0 .012%，誤差非常小。因此該DBM 能夠有效模擬爆轟現象。

圖7 爆轟波附近的物理量分布Fig. 7 Profiles of physical quantities around the detonation wave

3 結 論

基于動理學理論，本文提出了適用于超聲速可壓縮化學反應流的二維簡化DBM。該模型具有下列特點：

1）構建并使用了僅有9 個離散速度的模型，即離散速度分為3 組，每組含有3 個速度，各組速度的方向均勻分布、大小獨立可調。由于離散速度較少，模型具有較高的運算效率。

2）為了描述分子轉動和振動對應的額外自由度，引入了三組獨立可調的參數用于描述額外自由度部分的內能。由此，該DBM 具備了模擬比熱比可調的反應流系統的功能。

3）平衡態離散速度分布函數與化學反應項各自都滿足9 個獨立的動理學矩關系，并可以通過矩陣求逆的方式獲得其數學表達式。該方法具有物理精確、運算高效的特點。

4）通過Chapman-Enskog 多尺度分析可以證明，該DBM 除了在連續性極限條件下可以恢復包含化學反應的Euler 方程組之外，還具有描述一定熱力學非平衡行為的功能。

5）使用形式統一的離散玻爾茲曼方程，方法簡單、易于編程。同時，離散玻爾茲曼演化方程具有時空局域性特點，易于實現并行化，方便在高性能計算集群運行。

本文開展了3 個數值算例，驗證了該DBM 的可靠性。首先，通過均勻反應系統驗證模型能夠將化學反應與物理場自然耦合，適用于吸熱或放熱反應，并且反應流系統的比熱比可調；其次，通過激波管驗證模型適用于可壓縮流體系統，并且能夠同時捕捉稀疏波、物質界面和沖擊波的時空演化；最后，通過爆轟波算例開展了網格無關性驗證，并驗證模型適用于超聲速可壓縮化學反應流。

已有的二維DBM 至少包含16 個離散速度、動理學矩關系和演化方程，等價于N-S（或者Burnett）方程組和包含更多非平衡效應的粗粒化模型[24,30,36]。與之相比，簡化DBM 的優勢為離散速度數量較少、計算效率較高；局限性為物理精度較低，適用于熱力學非平衡程度較弱的物理系統。