高速可壓反應(yīng)流的二維簡(jiǎn)化離散玻爾茲曼模型

林傳棟

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

0 引 言

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

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

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

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

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

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

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

目前，DBM 在化學(xué)反應(yīng)流方面已經(jīng)得到初步應(yīng)用。2013 年，DBM 已被初步發(fā)展到爆轟領(lǐng)域，并被用于模擬分析爆轟波附近的非平衡效應(yīng)[22]。2014 年，林傳棟等[23]構(gòu)建并使用極坐標(biāo)系DBM 模擬研究了內(nèi)爆和外爆過程的非平衡強(qiáng)度，并從動(dòng)理學(xué)角度分析了粒子速度分布函數(shù)的主要特征。2015 年，許愛國等[24]構(gòu)建并使用多松弛DBM 模擬研究了黏性和熱傳導(dǎo)在爆轟系統(tǒng)中對(duì)局域和全域非平衡效應(yīng)的影響。2016 年，張玉東等[26]使用DBM 模擬研究了負(fù)溫度系數(shù)對(duì)爆轟系統(tǒng)的影響機(jī)理。之后，林傳棟等構(gòu)建了多組分DBM 用于模擬預(yù)混、非預(yù)混和部分預(yù)混反應(yīng)流系統(tǒng)[27]，又使用高效、準(zhǔn)確的多松弛DBM研究了不穩(wěn)定反應(yīng)熱流演化過程的非平衡行為[30]。

為了在連續(xù)性極限條件下恢復(fù)N-S 方程組或者更高階的流體力學(xué)方程組，并且為了包含更準(zhǔn)確和詳實(shí)的熱力學(xué)非平衡信息，要求DBM 包含足夠多的動(dòng)理學(xué)矩關(guān)系?，F(xiàn)有的DBM 模型至少能夠恢復(fù)N-S方程組，與之對(duì)應(yīng)，至少包含了7 個(gè)張量形式的動(dòng)理學(xué)矩關(guān)系。例如，對(duì)于二維DBM，當(dāng)包含額外自由度時(shí)，需要構(gòu)建16 速度的離散速度模型[30]，與之對(duì)應(yīng)，需要包含16 組形式統(tǒng)一的DBM 演化方程。然而，在實(shí)際的物理系統(tǒng)中，有時(shí)（或有些區(qū)域）物理場(chǎng)蘊(yùn)含的非平衡效應(yīng)可能較弱，不需要特精確和太費(fèi)時(shí)的模擬工具，而需要使用較為簡(jiǎn)潔和高效的模型。

為此，本文首次提出并構(gòu)建了簡(jiǎn)化的二維DBM。該模型只包含9 個(gè)獨(dú)立的動(dòng)理學(xué)矩關(guān)系，并且在連續(xù)性極限條件下能夠恢復(fù)帶有化學(xué)反應(yīng)的Euler 方程組。該DBM 包含9 個(gè)離散速度，基于9 組形式統(tǒng)一的DBM 演化方程。相對(duì)于之前構(gòu)建的16 速度模型[30]或者24 速度模型[24]，該DBM 具有更高的運(yùn)算效率。本項(xiàng)目的開展有助于推動(dòng)動(dòng)理學(xué)方法在化學(xué)反應(yīng)流的發(fā)展應(yīng)用，并進(jìn)一步促進(jìn)其在復(fù)雜多尺度物理系統(tǒng)的模擬研究。

1 二維簡(jiǎn)化DBM

1.1 離散玻爾茲曼方程

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

1.2 動(dòng)理學(xué)矩關(guān)系

1.3 化學(xué)反應(yīng)進(jìn)程

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

1.4 流體力學(xué)模型和非平衡物理量

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

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

另外，在已有的二維DBM 中，平衡態(tài)分布函數(shù)和反應(yīng)項(xiàng)都至少各自滿足16 個(gè)[30]（或24 個(gè)[24]）獨(dú)立的矩關(guān)系，DBM 等價(jià)于帶有化學(xué)反應(yīng)的N-S 方程組和包含較多非平衡效應(yīng)的粗?；Ｐ?。尤其對(duì)于文獻(xiàn)[36]，平衡態(tài)分布函數(shù)滿足26 個(gè)獨(dú)立的矩關(guān)系，該DBM 等價(jià)于Burnett 方程組和包含更多非平衡效應(yīng)的粗?；Ｐ?。在本文構(gòu)建的簡(jiǎn)化模型中，平衡態(tài)分布函數(shù)和反應(yīng)項(xiàng)都僅滿足9 個(gè)獨(dú)立的矩關(guān)系。相比于已有的DBM，該DBM 的局限性為：滿足的動(dòng)理學(xué)矩關(guān)系較少、物理精度較低，僅適用于熱力學(xué)非平衡程度較弱的物理系統(tǒng)。

1.5 數(shù)值離散格式

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

1.5.1 Runge-Kutta 格式

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

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

1）程序的開頭是初始化，根據(jù)物理場(chǎng)的初始條件，輸入?yún)?shù)；

2）根據(jù)給定的時(shí)間步數(shù)，即程序循環(huán)周期，判斷是否終止運(yùn)行；

3）根據(jù)給定的時(shí)間間隔，判斷是否輸出數(shù)值結(jié)果；

需要說明的是，上述NND 格式具有二階空間精度，是二階中心格式和二階迎風(fēng)格式的結(jié)合，是具有負(fù)系數(shù)的四階耗散項(xiàng)。因此，該NND 格式可以抑制奇偶失聯(lián)震蕩，能夠有效捕捉?jīng)_擊波和爆轟波[38]。

2 數(shù)值驗(yàn)證

本部分通過3 個(gè)典型算例對(duì)DBM 進(jìn)行數(shù)值驗(yàn)證。首先，通過均勻化學(xué)反應(yīng)系統(tǒng)驗(yàn)證本模型能夠?qū)⒒瘜W(xué)反應(yīng)與物理場(chǎng)自然耦合，模型既適用于吸熱反應(yīng)，也適用于放熱反應(yīng)，并且反應(yīng)流系統(tǒng)的比熱比可調(diào)；其次，通過Sod 激波管驗(yàn)證本模型適用于可壓縮流體系統(tǒng)，并且能夠同時(shí)捕捉稀疏波、物質(zhì)界面和沖擊波的時(shí)空演化；最后，通過爆轟波算例開展了網(wǎng)格無關(guān)性驗(yàn)證，并驗(yàn)證本模型適用于超聲速可壓縮化學(xué)反應(yīng)流。

2.1 均勻化學(xué)反應(yīng)系統(tǒng)

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

圖3 展示了化學(xué)反應(yīng)后不同化學(xué)能對(duì)應(yīng)的系統(tǒng)溫度，其關(guān)系見式（37），此處設(shè)置比熱比 γ=1.4。符號(hào)代表DBM 模擬結(jié)果，實(shí)線代表理論解：

可以看出，DBM 的模擬結(jié)果與理論解完全吻合。而且，圖3 中化學(xué)能的數(shù)值涵蓋負(fù)值到正值的范圍。其中負(fù)值區(qū)域代表化學(xué)反應(yīng)吸熱，正值區(qū)域代表化學(xué)反應(yīng)放熱。也就是說，該DBM 既適用于化學(xué)吸熱系統(tǒng)，也適用于化學(xué)放熱系統(tǒng)。

圖3 化學(xué)反應(yīng)后不同化學(xué)能對(duì)應(yīng)的系統(tǒng)溫度Fig. 3 Temperature after chemical reaction under various chemical energies

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

圖4 化學(xué)反應(yīng)后不同比熱比對(duì)應(yīng)的系統(tǒng)溫度Fig. 4 Temperature after chemical reaction under various specific heat ratios

2.2 激波管

圖5 展示了在t= 0.02 時(shí)刻激波管中的物理量分布，圖5（a～d）依次是密度、溫度、壓強(qiáng)和水平速度。符號(hào)代表DBM 模擬結(jié)果，實(shí)線代表理論解?？梢钥闯觯涸谧钭髠?cè)的是稀疏波（向左傳播），在中間的是物質(zhì)分界面，在最右側(cè)的是沖擊波（向右傳播）。其中，稀疏波在水平方向的跨度最大，各物理量都有清晰的空間梯度；在物質(zhì)界面兩側(cè)，密度和溫度具有明顯空間梯度，而壓強(qiáng)和速度均勻分布；在沖擊波波陣面，各物理量梯度最大，具有劇烈的空間變化。在以上演化過程中，DBM 與理論解具有很高的匹配度。這充分說明，DBM 能夠很好地捕捉稀疏波、物質(zhì)界面和沖擊波的空間分布。

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

2.3 爆轟波

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

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

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

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

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

圖7 展示了在空間步長Δx=1×10?4、 時(shí)間步長Δt=2×10?6下爆轟波附近的模擬結(jié)果。圖7（a～d）依次對(duì)應(yīng)密度、溫度、壓強(qiáng)和水平速度的空間分布，其中符號(hào)代表DBM 模擬結(jié)果，實(shí)線代表Zeldovich-Neumann-Doering（ZND）理論解[1]?？梢钥闯?，DBM模擬結(jié)果與ZND 理論解具有非常好的一致性。尤其是當(dāng)化學(xué)反應(yīng)結(jié)束后，在爆轟波后側(cè)形成穩(wěn)定的“平臺(tái)區(qū)域”，此處的DBM 模擬結(jié)果為 ρ=1.48065、T=2.06216、p=3.05334、ux=?1.69974。對(duì)比ZND 理論解，相對(duì)誤差依次為： 0 .015%、 0 .047%、0 .032%、0 .012%，誤差非常小。因此該DBM 能夠有效模擬爆轟現(xiàn)象。

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

3 結(jié) 論

基于動(dòng)理學(xué)理論，本文提出了適用于超聲速可壓縮化學(xué)反應(yīng)流的二維簡(jiǎn)化DBM。該模型具有下列特點(diǎn)：

1）構(gòu)建并使用了僅有9 個(gè)離散速度的模型，即離散速度分為3 組，每組含有3 個(gè)速度，各組速度的方向均勻分布、大小獨(dú)立可調(diào)。由于離散速度較少，模型具有較高的運(yùn)算效率。

2）為了描述分子轉(zhuǎn)動(dòng)和振動(dòng)對(duì)應(yīng)的額外自由度，引入了三組獨(dú)立可調(diào)的參數(shù)用于描述額外自由度部分的內(nèi)能。由此，該DBM 具備了模擬比熱比可調(diào)的反應(yīng)流系統(tǒng)的功能。

3）平衡態(tài)離散速度分布函數(shù)與化學(xué)反應(yīng)項(xiàng)各自都滿足9 個(gè)獨(dú)立的動(dòng)理學(xué)矩關(guān)系，并可以通過矩陣求逆的方式獲得其數(shù)學(xué)表達(dá)式。該方法具有物理精確、運(yùn)算高效的特點(diǎn)。

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

5）使用形式統(tǒng)一的離散玻爾茲曼方程，方法簡(jiǎn)單、易于編程。同時(shí)，離散玻爾茲曼演化方程具有時(shí)空局域性特點(diǎn)，易于實(shí)現(xiàn)并行化，方便在高性能計(jì)算集群運(yùn)行。

本文開展了3 個(gè)數(shù)值算例，驗(yàn)證了該DBM 的可靠性。首先，通過均勻反應(yīng)系統(tǒng)驗(yàn)證模型能夠?qū)⒒瘜W(xué)反應(yīng)與物理場(chǎng)自然耦合，適用于吸熱或放熱反應(yīng)，并且反應(yīng)流系統(tǒng)的比熱比可調(diào)；其次，通過激波管驗(yàn)證模型適用于可壓縮流體系統(tǒng)，并且能夠同時(shí)捕捉稀疏波、物質(zhì)界面和沖擊波的時(shí)空演化；最后，通過爆轟波算例開展了網(wǎng)格無關(guān)性驗(yàn)證，并驗(yàn)證模型適用于超聲速可壓縮化學(xué)反應(yīng)流。

已有的二維DBM 至少包含16 個(gè)離散速度、動(dòng)理學(xué)矩關(guān)系和演化方程，等價(jià)于N-S（或者Burnett）方程組和包含更多非平衡效應(yīng)的粗粒化模型[24,30,36]。與之相比，簡(jiǎn)化DBM 的優(yōu)勢(shì)為離散速度數(shù)量較少、計(jì)算效率較高；局限性為物理精度較低，適用于熱力學(xué)非平衡程度較弱的物理系統(tǒng)。