許愛國，單奕銘, 陳鋒，甘延標(biāo)，林傳棟

1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所 計(jì)算物理實(shí)驗(yàn)室，北京 100088 2.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081 3.北京大學(xué) 應(yīng)用物理與技術(shù)研究中心 高能量密度物理數(shù)值模擬教育部重點(diǎn)實(shí)驗(yàn)室，北京 100871 4.山東交通學(xué)院 航空學(xué)院，濟(jì)南 250357 5.北華航天工業(yè)學(xué)院 河北省跨氣水介質(zhì)飛行器重點(diǎn)實(shí)驗(yàn)室，廊坊 065000 6.中山大學(xué) 中法核工程與技術(shù)學(xué)院，珠海 519082

隨著科技進(jìn)步，人類飛行器的速度越來越高?！案咚贇饬髦械狞c(diǎn)火與持續(xù)燃燒成為空天推進(jìn)中長(zhǎng)久的關(guān)鍵問題”[1]。這些關(guān)鍵問題，往往涉及極端條件。極端條件，是相對(duì)于常溫常壓常速等普通條件而言的，它至少在某一方面嚴(yán)重偏離普通條件。隨著科技發(fā)展，超聲速、高瞬變已經(jīng)成為航空領(lǐng)域司空見慣的條件和工況[2-3]。美國在20世紀(jì) 末提出了Hyper-X研究計(jì)劃，通過研究一體化超燃沖壓發(fā)動(dòng)機(jī)技術(shù)從而實(shí)現(xiàn)馬赫數(shù)7～15的飛行[4]；俄羅斯也一直在進(jìn)行高超聲速飛行器的研究工作，比如“冷”計(jì)劃(馬赫數(shù)達(dá)到6.5，采用超燃沖壓發(fā)動(dòng)機(jī))、“鷹”計(jì)劃(馬赫數(shù)達(dá)到12～14)、“針”計(jì)劃(馬赫數(shù)達(dá)到6～14，為空天飛機(jī)樣機(jī))。飛行器的高速運(yùn)動(dòng)導(dǎo)致流體與壁面間存在著極強(qiáng)的相互作用，大量動(dòng)能轉(zhuǎn)變?yōu)闊崮?，高溫使得理想氣體模型失效，真實(shí)氣體效應(yīng)影響顯著，同時(shí)也會(huì)對(duì)飛行器表面產(chǎn)生燒蝕作用，使其防熱結(jié)構(gòu)遭到破壞，從而改變飛行器氣動(dòng)性能，影響其穩(wěn)定性和飛行安全性。并且在高溫下飛行器表面的防熱材料也會(huì)影響邊界層內(nèi)氣體化學(xué)反應(yīng)的催化效應(yīng)，不同壁面催化條件會(huì)對(duì)氣動(dòng)熱影響很大。

超聲速、高瞬變等極端條件引發(fā)一系列值得深入研究的可壓、非平衡復(fù)雜流動(dòng)行為；非線性效應(yīng)、離散效應(yīng)、強(qiáng)耦合效應(yīng)是其重要特征。這給傳統(tǒng)流體建模理論的合理性帶來嚴(yán)重挑戰(zhàn)；同時(shí)，這些行為中的很大一部分又發(fā)生在微觀分子動(dòng)力學(xué)模擬無法企及的時(shí)間和空間尺度。例如，高空環(huán)境中的飛行器周圍空氣密度低，稀薄大氣條件下連續(xù)介質(zhì)假設(shè)失效，并且在不同的飛行高度上高超聲速飛行器所處的環(huán)境狀態(tài)也是不同的，由高空到低空分別處于自由分子流、過渡流、滑移流和連續(xù)流狀態(tài)，這就需要對(duì)飛行進(jìn)行跨尺度模擬。氣體經(jīng)過飛行器頭部的強(qiáng)激波被極速加熱，在氣體稀薄效應(yīng)和高溫效應(yīng)的共同作用下，氣體分子的平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)、電子能被同時(shí)激發(fā)，化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間幾乎相當(dāng)，飛行器的外流場(chǎng)表現(xiàn)出明顯的非平衡性。而且，飛行器在飛行中，其周圍的氣體還存在層流到湍流的轉(zhuǎn)捩現(xiàn)象，在湍流邊界層的氣動(dòng)熱更加嚴(yán)重。

基于求解Navier-Stokes方程組的傳統(tǒng)計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)已經(jīng)在許多領(lǐng)域取得巨大的成功，但在航空、航天、微流控等領(lǐng)域也遇到了諸多新的瓶頸與挑戰(zhàn)。其原因分為2個(gè)方面：① 物理建模層面的問題；② 離散格式帶來的數(shù)值精度和穩(wěn)定性問題。物理模型合理和具備相應(yīng)功能是數(shù)值仿真研究的前提；物理建模層面的問題是無法通過數(shù)值精度的提高來解決的。鑒于一些燃燒相關(guān)的概念和技術(shù)(例如脈沖和旋轉(zhuǎn)爆轟、微尺度燃燒、納米推進(jìn)劑、部分預(yù)混和層流燃燒、等離子體燃燒、冷火焰等)均涉及非平衡問題[5-19]，也已有大量文獻(xiàn)討論離散格式等數(shù)值問題，同時(shí)鑒于作者團(tuán)隊(duì)的知識(shí)結(jié)構(gòu)和研究經(jīng)驗(yàn)主要在流體物理方面，本文從物理建模與復(fù)雜物理場(chǎng)分析角度，介紹離散玻爾茲曼建模方法(Discrete Boltzmann Modeling method, DBM)的研究進(jìn)展[20-28]。本文結(jié)構(gòu)如下：第1節(jié)回顧宏觀與微觀建模方法的特點(diǎn)；第2節(jié)介紹燃燒系統(tǒng)的DBM建模理論；第3節(jié)給出部分?jǐn)?shù)值模擬結(jié)果；最后給出結(jié)論與說明。

1 宏觀與微觀建模

人們先前對(duì)燃燒和爆轟問題的研究主要依賴于實(shí)驗(yàn)及少量的理論分析[29-30]，隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，燃燒和爆轟的數(shù)值模擬研究取得了巨大的成就[31-38]。數(shù)值模擬基于物理建模。燃燒系統(tǒng)的物理建模，有3個(gè)層次：微觀、介觀和宏觀尺度。微觀層次通常是指分子動(dòng)力學(xué)(Molecular Dynamics, MD)描述[39-41]，其中分子間作用勢(shì)的確立是模型構(gòu)建的關(guān)鍵，分子間作用勢(shì)有效半徑的選取是模擬過程中的關(guān)鍵，有效半徑的選取通常以滿足需求的最小值為最佳原則。微觀研究可以幫助建立反應(yīng)速率方程，提供完備信息，但由于運(yùn)算量太大，對(duì)內(nèi)存要求太高，因而適用的時(shí)空尺度非常有限。宏觀尺度通常是指基于Euler或Navier-Stokes方程外加一個(gè)化學(xué)反應(yīng)唯象模型的描述，其模擬手段以傳統(tǒng)CFD方法為主[5,42-43]。近年來格子玻爾茲曼方法(Lattice Boltzmann Method, LBM)正在引起廣泛興趣[44-53]，成為求解燃燒系統(tǒng)宏觀流體力學(xué)方程組的另一手段。在宏觀層次上人們關(guān)注的主要是系統(tǒng)的近平衡宏觀行為，即密度、溫度、流速、壓強(qiáng)、應(yīng)力、熱流等物理量的演化，控制方程是代表質(zhì)量、動(dòng)量和能量守恒的流體力學(xué)方程組，本構(gòu)關(guān)系往往是根據(jù)經(jīng)驗(yàn)或唯象理論給出的。

2 介觀動(dòng)理學(xué)建模理論

從動(dòng)理學(xué)理論和Chapman-Enskog多尺度分析角度看，如果系統(tǒng)時(shí)刻均處于局域熱動(dòng)平衡態(tài)，即在Chapman-Enskog多尺度分析時(shí)只保留Knudsen數(shù)(Kn)的零次方項(xiàng)，則Boltzmann方程對(duì)應(yīng)的宏觀流體方程組就是Euler方程組；如果系統(tǒng)時(shí)刻處于熱動(dòng)平衡態(tài)附近，二階及以上非平衡效應(yīng)很弱以至于可以忽略，只需考慮一階非平衡效應(yīng)，那么此時(shí) Boltzmann方程對(duì)應(yīng)的宏觀流體方程組是Navier-Stokes方程組。也就是說，Euler方程實(shí)際上假設(shè)系統(tǒng)始終處于熱力學(xué)平衡態(tài)，Navier-Stokes方程只包含了一階的非平衡效應(yīng)。而包括燃燒、爆轟在內(nèi)的復(fù)雜流體系統(tǒng)內(nèi)部往往存在大量的中間尺度的結(jié)構(gòu)和動(dòng)理學(xué)模式。結(jié)構(gòu)小到一定程度會(huì)引發(fā)“離散效應(yīng)”(相對(duì)于所關(guān)注的結(jié)構(gòu)尺度，平均分子間距不再是可以忽略的小量，不能再使用連續(xù)性假設(shè))，模式快到一定程度會(huì)引發(fā)“熱力學(xué)非平衡效應(yīng)”(相對(duì)于所關(guān)注的模式的時(shí)間尺度，熱力學(xué)弛豫時(shí)間不再是可以忽略的小量，不能再假設(shè)在模式演化的每一步系統(tǒng)都已回到熱力學(xué)平衡態(tài))。因而，系統(tǒng)內(nèi)的這些小結(jié)構(gòu)、快模式行為向Navier-Stokes方程等傳統(tǒng)流體力學(xué)理論的物理合理性提出了挑戰(zhàn)。同時(shí)，人們關(guān)注的行為特征往往又發(fā)生在微觀分子動(dòng)力學(xué)模擬無法企及的時(shí)空尺度上。這些“介尺度”的結(jié)構(gòu)和行為因?yàn)槟Ｐ秃头椒ǖ娜狈?不成熟)而研究薄弱。這些薄弱不僅在一定程度上阻斷了對(duì)微觀到宏觀之間關(guān)聯(lián)的認(rèn)識(shí)，也在一定程度上阻礙了本構(gòu)模型的機(jī)理化、科學(xué)化；同時(shí)，微流控等技術(shù)的快速發(fā)展在提醒人們，這些(相對(duì)于宏觀)特征更加豐富但以前知之甚少的“介尺度”非平衡、不連續(xù)行為，往往意味著大量待開發(fā)的物理功能。離散玻爾茲曼建模方法就是在這個(gè)背景下產(chǎn)生的理論模型構(gòu)建方法。它是統(tǒng)計(jì)物理學(xué)相空間描述方法在離散玻爾茲曼方程形式下的進(jìn)一步發(fā)展[54-59]，其思想起源于許愛國等于2012年發(fā)表的一篇研究綜述[54]；在發(fā)展過程中又受到了形態(tài)相空間描述方法的啟發(fā)[59-60]。DBM的研究思路是：將復(fù)雜問題進(jìn)行分解，根據(jù)研究需求，選取一個(gè)視角，研究系統(tǒng)的一組動(dòng)理學(xué)性質(zhì)，因而要求描述這組性質(zhì)的動(dòng)理學(xué)矩在模型簡(jiǎn)化過程中保值；以該組動(dòng)理學(xué)矩的獨(dú)立分量為基，構(gòu)建相空間，使用該相空間和其子空間來描述系統(tǒng)的狀態(tài)和行為；研究視角和建模精度隨著研究推進(jìn)和實(shí)際需求而調(diào)整。一個(gè)DBM模型的建立，需要經(jīng)歷3個(gè)步驟。第1步，引入一個(gè)形式上的局域目標(biāo)分布函數(shù)，將原來的碰撞項(xiàng)寫成一個(gè)線性化碰撞算符的形式；要求是，所關(guān)心的物理特征量使用簡(jiǎn)化前和簡(jiǎn)化后的模型計(jì)算，所得的結(jié)果必須一致。第2步，借助離散速度，將原本連續(xù)、積分形式的動(dòng)理學(xué)矩轉(zhuǎn)化為求和進(jìn)行計(jì)算；要求是，所關(guān)心的動(dòng)理學(xué)矩轉(zhuǎn)換為求和進(jìn)行計(jì)算后，得到的結(jié)果必須相同，即

(1)

式中：ψ(v)=[1,v,vv,…]對(duì)應(yīng)研究中關(guān)心的、建模過程中要保值的動(dòng)理學(xué)矩，v為分子速度；fi=f(vi)，vi為第i個(gè)離散速度。第3步，是DBM建模的目的和核心，給出非平衡狀態(tài)和行為描述的具體方案。DBM建模主要針對(duì)的是宏觀連續(xù)建模失效或物理功能不足、而微觀分子動(dòng)力學(xué)模擬又因?yàn)檫m用尺度問題無能為力的“介尺度”“兩難”情形。DBM提供的物理信息量介于宏觀連續(xù)描述和微觀分子動(dòng)力學(xué)之間。相對(duì)于宏觀描述，DBM從一個(gè)更寬的視角來觀測(cè)系統(tǒng)；DBM中非守恒矩描述的必要性和收益均隨著非平衡程度的增高而增加。其與格子玻爾茲曼方法[44]、格子氣方法[61]的關(guān)系如圖1所示。早期的格子氣方法，其本身就孕育了(至少)2個(gè)發(fā)展方向：① 統(tǒng)計(jì)物理學(xué)領(lǐng)域的粗?；７椒?；② 計(jì)算數(shù)學(xué)領(lǐng)域的方程解法。1952年，李政道和楊振寧發(fā)表在《Physical Review》的《Statistical theory of equations of state and phase transitions.II.Lattice gas and Ising model》便是統(tǒng)計(jì)物理學(xué)領(lǐng)域構(gòu)建和使用格子氣粗?；７椒ǖ膶?shí)例之一[62]。通常認(rèn)為，1988年開始出現(xiàn)LBM方法的雛形[61]。在1988—2012年期間，LBM方程解法和LBM建模方法在相互交織、相互啟發(fā)中發(fā)展。LBM建模方法呈現(xiàn)出來的主要功能還是恢復(fù)宏觀流體模型；與LBM方程解法相比，只是其在遵守物理規(guī)則方面要求更加嚴(yán)格，例如不允許使用非統(tǒng)計(jì)物理學(xué)意義下的“玻爾茲曼方程”、“矩關(guān)系”等，并未表現(xiàn)出與求解流體方程(組)功能的顯著差異。于是，LBM基本上成了LBM方程解法的代名詞。2012年，LBM方法與統(tǒng)計(jì)物理學(xué)非平衡行為基本描述方法相遇并碰撞出火花，LBM方法被注入統(tǒng)計(jì)物理學(xué)使用分布函數(shù)非守恒矩來描述非平衡狀態(tài)和行為的思路，或者說統(tǒng)計(jì)物理學(xué)非平衡行為描述方法在離散玻爾茲曼方程形式下獲得進(jìn)一步發(fā)展[54]。隨后，LBM建模方法與LBM方程解法的功能差異開始逐漸變得清晰。其中，起到關(guān)鍵促進(jìn)作用的一步是相空間描述方法在離散玻爾茲曼方程形式下的進(jìn)一步發(fā)展，這在后面還要提到。鑒于經(jīng)常被誤解為L(zhǎng)BM方程解法，LBM建模方法逐漸改稱為L(zhǎng)B動(dòng)理學(xué)建模(LB Kinetic Modeling, LBKM)、DBM。圖1中DBM建模與分析方法軸上2012以后的點(diǎn)對(duì)應(yīng)燃燒系統(tǒng)DBM建模文章發(fā)表的年份。需要指出的是，統(tǒng)計(jì)物理學(xué)領(lǐng)域作為粗粒化建模的格子氣方法，不會(huì)因?yàn)镈BM或其他建模方法的出現(xiàn)而消失，它的思路永遠(yuǎn)相對(duì)獨(dú)立地存在和發(fā)展，在需要它的地方發(fā)揮作用。

圖1 DBM、形態(tài)學(xué)描述方法和LBM的發(fā)展歷程

模擬燃燒系統(tǒng)的DBM，其演化方程可統(tǒng)一寫為[56]

(2)

(3)

式中：λ為化學(xué)反應(yīng)進(jìn)程參數(shù)；F(λ)表示反應(yīng)速率函數(shù)。這樣，便得到了Ci的表達(dá)式為

(4)

(5)

其中：ρ為密度；u為流速；Q為單位質(zhì)量反應(yīng)物完全反應(yīng)后釋放的熱量，在本文中稱之為爆熱；e為內(nèi)能密度。

(6)

(7)

圖2 非平衡特征量張開的相空間

盡管爆轟研究已有100多年的歷史, 但時(shí)至今日，它仍然是國際熱點(diǎn)研究問題之一[11,13-14,16,66-69]。到目前為止，幾乎所有獲得廣泛應(yīng)用的化學(xué)反應(yīng)模型均是唯象的或半唯象的。例如, Arrhenius反應(yīng)率、森林火災(zāi)模型、兩步模型、Cochran反應(yīng)率模型、Lee-Tarver模型等。在實(shí)際應(yīng)用過程中需要根據(jù)具體問題選擇合適的化學(xué)反應(yīng)模型。具體工況不同，點(diǎn)火的具體處理方式也就可能不同：化學(xué)反應(yīng)的啟動(dòng)可能是由沖擊引起的，也可能是烤燃、電擊或摩擦等引起的，但其底層物理原因均是由于溫度超過了化學(xué)反應(yīng)的啟動(dòng)溫度，反應(yīng)物分子內(nèi)的原子或原子團(tuán)通過熱運(yùn)動(dòng)掙脫它們之間的化學(xué)勢(shì)引起的。在已發(fā)表的DBM文獻(xiàn)中，使用的化學(xué)反應(yīng)啟動(dòng)判據(jù)均是溫度判據(jù)。同時(shí)需要說明的是，DBM既可以用于研究爆轟，也可以用于研究火焰。但一般意義下的火焰引發(fā)的熱力學(xué)非平衡強(qiáng)度往往弱于爆轟，特別是強(qiáng)爆轟。非平衡效應(yīng)越弱，則已有的傳統(tǒng)流體建模就越能勝任；非平衡效應(yīng)越強(qiáng)，則傳統(tǒng)流體建模在物理功能方面就呈現(xiàn)出更多不足，DBM的必要性和優(yōu)勢(shì)也就越明顯。

3 部分研究結(jié)果

近年來，DBM在燃燒(含爆轟)建模與模擬方面取得了一系列的進(jìn)展，為非平衡燃燒(含爆轟)研究帶來一系列新的物理認(rèn)知和思考。

2014年，林傳棟等[70]基于極坐標(biāo)系構(gòu)建了一個(gè)用于模擬燃燒的DBM，模擬了典型的內(nèi)爆和外爆過程，觀測(cè)并分析了化學(xué)反應(yīng)熱、物質(zhì)輸運(yùn)、熱傳導(dǎo)與幾何(匯聚或發(fā)散)效應(yīng)的合作與競(jìng)爭(zhēng)。在外爆過程中觀測(cè)到了與一維爆轟情形不同的、幾何效應(yīng)導(dǎo)致的現(xiàn)象，例如熄爆、雙向爆轟、穩(wěn)定爆轟等?；瘜W(xué)反應(yīng)釋放的熱量使得系統(tǒng)局域溫度升高, 而熱傳導(dǎo)和幾何發(fā)散效應(yīng)使得系統(tǒng)局域溫度降低。所以，預(yù)先起爆區(qū)域的大小，影響著釋放的熱量是否足以克服幾何發(fā)散效應(yīng)引發(fā)的溫度降低。如果之前的化學(xué)反應(yīng)熱足夠多, 則化學(xué)反應(yīng)得以維持；如果熱傳導(dǎo)和幾何發(fā)散效應(yīng)占優(yōu)勢(shì), 則會(huì)出現(xiàn)熄爆現(xiàn)象。隨著爆轟波向外傳播, 幾何發(fā)散效應(yīng)逐漸減弱, 發(fā)生熄爆的可能性降低。發(fā)現(xiàn)在高度對(duì)稱的系統(tǒng)中，球心處系統(tǒng)始終處于熱力學(xué)平衡態(tài)。2018年，許愛國等提出一個(gè)單弛豫時(shí)間球坐標(biāo)DBM模型[71]。在內(nèi)爆與外爆過程中，幾何匯聚與發(fā)散效應(yīng)起到一個(gè)“外場(chǎng)力”的作用，更加細(xì)節(jié)的討論參見文獻(xiàn)[71]。

上述DBM均為單流體模型，即忽略介質(zhì)成分差異，只關(guān)注其平均行為，產(chǎn)物和反應(yīng)物所占份額用一個(gè)反應(yīng)進(jìn)程參數(shù)λ來描述。為了更加細(xì)致地描述反應(yīng)系統(tǒng)，例如可以同時(shí)觀測(cè)反應(yīng)物和產(chǎn)物的不同流速和溫度，2016年林傳棟等[72]提出一個(gè)二流體燃燒DBM模型，在該模型中所有的反應(yīng)物視為同一種介質(zhì)(成分)，所有的產(chǎn)物視為同一種介質(zhì)(成分)，反應(yīng)物和產(chǎn)物分別使用2個(gè)不同的分布函數(shù)來描述，用2個(gè)相耦合的離散玻爾茲曼方程來描述反應(yīng)物和產(chǎn)物的演化過程。該模型可用于模擬亞聲速和超聲速燃燒系統(tǒng)，比熱比可調(diào)。相關(guān)宏觀流體模型(例如帶有化學(xué)反應(yīng)的Navier-Stokes方程，F(xiàn)ick第一、第二定律，Stefan-Maxwell擴(kuò)散方程等)均是該模型的特例，借助該模型可以很方便地觀測(cè)和研究與宏觀行為相伴隨的各種熱力學(xué)非平衡效應(yīng)。為了更加細(xì)致地描述化學(xué)反應(yīng)系統(tǒng)，體現(xiàn)反應(yīng)物內(nèi)不同介質(zhì)(成分)、產(chǎn)物內(nèi)不同介質(zhì)(成分)之間的差異，進(jìn)一步提出一個(gè)多流體DBM模型。作為應(yīng)用實(shí)例，模擬研究了丙烷的氧化過程：C3H8+5O2→3CO2+4H2O。作為模型驗(yàn)證，除了密度、流速、溫度、壓強(qiáng)與理論解進(jìn)行比對(duì)，定壓燃燒的火焰溫度(T)與實(shí)驗(yàn)結(jié)果符合較好[73]，如圖3所示。因?yàn)椴煌姆磻?yīng)物得以分別描述，所以該模型既可以描述預(yù)混燃燒，又可以描述非預(yù)混燃燒。

2.7 準(zhǔn)確度考察（加樣回收率試驗(yàn)） 取錐形瓶6個(gè)，分別精密加入新配制的大黃素-8-O-β-D-葡萄糖苷對(duì)照品溶液（66 μg∕mL）6.36 mL、大黃素甲醚-8-O-β-D-葡萄糖苷對(duì)照品溶液（48 μg∕mL）3.75 mL，減壓回收溶劑至干，再精密稱取已知含量供試品（6號(hào)供試品粉末，大黃素-8-O-β-D-葡萄糖苷含量為4.2 mg∕g，大黃素甲醚-8-O-β-D-葡萄糖含量為1.8 mg∕g），各稱取約0.10 g，分別精密稱定，按“2.3”項(xiàng)下制備供試品溶液。進(jìn)樣測(cè)定，按下式計(jì)算回收率。

圖3 定壓燃燒火焰溫度分布[73]

關(guān)于含外場(chǎng)力情形的燃燒DBM建?？蓞⒖嘉墨I(xiàn)[74]。如果令式(2)中的Force Term=-Fi，則在分布函數(shù)偏離平衡部分(f-feq)對(duì)外場(chǎng)力效應(yīng)影響較小，且化學(xué)反應(yīng)的時(shí)間尺度遠(yuǎn)大于熱力學(xué)弛豫時(shí)間的情形，在BGK(Bhatanger-Gross-Krook)模型框架下外場(chǎng)力和化學(xué)反應(yīng)引發(fā)的分布函數(shù)變化率可寫為

(8)

式(8)的含義為：在時(shí)間間隔τ內(nèi)，流速由u變化為u+aτ，溫度由T變?yōu)門+τT′。由外力和化學(xué)反應(yīng)而引起的能量的變化率為

E′=ρu·a+ρQλ′

(9)

根據(jù)E=(D+I)ρT/2+ρu·u/2和式(9)，可以得到溫度的變化率為

(10)

式中：D代表維度；I代表與分子旋轉(zhuǎn)和/或內(nèi)部振動(dòng)相對(duì)應(yīng)的額外自由度。反應(yīng)進(jìn)度參數(shù)λ定義為反應(yīng)產(chǎn)物與混合物的質(zhì)量比。在文獻(xiàn)[74]中，作為實(shí)例，化學(xué)反應(yīng)進(jìn)程由Cochran速率函數(shù)控制：

λ′=ω1pm(1-λ)+ω2pnλ(1-λ)

(11)

由式(11)可知，λ′依賴于壓強(qiáng)p=ρT，而ω1、ω2、m和n均為可調(diào)參數(shù)。

圖4 爆轟波附近的各物理量[75]

圖5 反應(yīng)物和產(chǎn)物在壓強(qiáng)峰值處的速度分布函數(shù)主要特征[75]

圖6 非平衡效應(yīng)與化學(xué)反應(yīng)放熱的關(guān)系[75]

為了估算爆轟波波峰的相對(duì)高度，定義了峰高H(q)=(qmax-qs)/(qvon-qs)，其中q指代各宏觀物理量，qmax為在爆轟波附近q的最大值，qs為CJ解而qvon是在馮·紐曼峰處的ZND解。除此之外，通過比較可以發(fā)現(xiàn)：① 考慮非平衡效應(yīng)的爆轟波波峰要低于不考慮非平衡效應(yīng)的結(jié)果；② 考慮非平衡效應(yīng)的爆轟波波陣面比不考慮非平衡效應(yīng)情形要寬；③ 考慮非平衡效應(yīng)的物理量的梯度要比不考慮非平衡效應(yīng)情形的小。

2016年，張玉東等[23]利用DBM研究了帶有爆轟的反應(yīng)流問題，從理論上推導(dǎo)了一套新的流體力學(xué)方程，方程中的應(yīng)力和熱流用文中定義的2個(gè)非平衡量(無組織動(dòng)量流和無組織能量流)進(jìn)行代替?；谒岢龅膭?dòng)理學(xué)模型，將2個(gè)非平衡量和熵產(chǎn)生速率之間建立起關(guān)系，并研究了負(fù)溫度系數(shù)對(duì)爆轟演化過程的影響。使用的化學(xué)反應(yīng)率模型為

(12)

其中:

(13)

(14)

(15)

式中：Tth為起爆溫度；k為化學(xué)反應(yīng)速率系數(shù)；h1和h2分別代表了k的峰值和谷值；T1和T2分別是h1和h2對(duì)應(yīng)的溫度。研究了圖7(a)～圖7(d)所示4種反應(yīng)速率系數(shù)下的爆轟現(xiàn)象。通過比較可以發(fā)現(xiàn)，對(duì)于4種反應(yīng)速率特性的爆轟情形，在遠(yuǎn)離爆轟波陣面附近以及化學(xué)反應(yīng)區(qū)，無組織動(dòng)量流與應(yīng)力、無組織能量流與熱流均有一些差異，而這些區(qū)域也正是非平衡特征最顯著的地方。當(dāng)?shù)氐撵禺a(chǎn)生有3個(gè)來源：化學(xué)反應(yīng)、無組織動(dòng)量流以及無組織能量流。對(duì)于系統(tǒng)內(nèi)的全局熵產(chǎn)生，化學(xué)反應(yīng)所占的比例遠(yuǎn)大于另外2個(gè)方面，無組織動(dòng)量流導(dǎo)致的熵產(chǎn)生率大于無組織能量流產(chǎn)生的熵產(chǎn)生率。負(fù)溫度系數(shù)對(duì)動(dòng)力學(xué)量的作用是降低馮·紐曼峰的高度(減小密度、壓強(qiáng)和速度)，加寬反應(yīng)區(qū)，抑制化學(xué)反應(yīng)進(jìn)程。

圖7 5種反應(yīng)速率系數(shù)隨溫度變化特征圖[23]

2019年，張玉東等[76]提出了一個(gè)用于模擬爆轟的一維離散玻爾茲曼模型。通過對(duì)Sod激波管問題、Colella爆炸波問題和一維自持穩(wěn)定爆轟傳播問題進(jìn)行模擬，并與相應(yīng)的解析解進(jìn)行比較，數(shù)值驗(yàn)證了該模型的有效性。在該工作中研究了由于負(fù)溫度系數(shù)而導(dǎo)致的一種反常爆轟現(xiàn)象，此時(shí)的化學(xué)反應(yīng)速率系數(shù)k的表達(dá)式為

k(T)=a+b[T3/3-(T1+T2)T2/2+T1T2T]

(16)

式中：系數(shù)a和b由式(14)和式(15)給出。在本工作中，h1=2 000,h2=10，T1=1.14，T2=1.45，起爆溫度Tth=1.1。化學(xué)反應(yīng)速率系數(shù)與溫度之間的關(guān)系如圖7(e)所示，從圖中可以清楚地看到，存在一個(gè)負(fù)溫度系數(shù)區(qū)間(Ti，Tj)，反應(yīng)速率系數(shù)隨著溫度的升高而減小。

圖8[76]為恒定的化學(xué)反應(yīng)速率系數(shù)與負(fù)溫度系數(shù)2種情形下的爆轟波模擬結(jié)果。2個(gè)結(jié)果之間除了化學(xué)反應(yīng)速率和網(wǎng)格數(shù)之外其他參數(shù)均一致，從圖中可以明顯發(fā)現(xiàn)在負(fù)溫度系數(shù)條件下出現(xiàn)了反常的爆轟現(xiàn)象。對(duì)于反常爆轟來說，其并沒有一個(gè)恒定的波速，并且波形是隨時(shí)間周期性變化的。

圖8 正常爆轟波與反常爆轟波的比較[76]

為了方便討論，根據(jù)溫度將化學(xué)反應(yīng)大致分為3個(gè)階段，如圖7(e)所示，這3個(gè)階段分別記為S1、S2和S3，其中第1階段S1處于低溫區(qū)域，但由于負(fù)溫度系數(shù)的緣故導(dǎo)致其具有較快的反應(yīng)速率，第2階段S2在特定溫度范圍內(nèi)具有較慢的反應(yīng)速率，而第3階段S3在較高溫度下具有較快的反應(yīng)速率并且反應(yīng)速率隨著溫度的升高而顯著的增加。通過研究發(fā)現(xiàn)，對(duì)于常規(guī)的爆轟現(xiàn)象，化學(xué)反應(yīng)主要發(fā)生在S1和S2階段中，但對(duì)于該反常爆轟，在S3階段中的某個(gè)時(shí)間會(huì)產(chǎn)生一個(gè)局部熱點(diǎn)，然后在舊的爆轟波波陣面后方出現(xiàn)具有更劇烈化學(xué)反應(yīng)的新爆轟波，這個(gè)新爆轟波擁有比前面爆轟波更快的波速，它緊追前波，最后2個(gè)爆轟波合并，之后爆轟波的波速開始減慢直到它達(dá)到CJ爆轟值。在此之后局部熱點(diǎn)再次出現(xiàn)并重復(fù)上述過程，如圖9[76]所示。

圖9 反常爆轟現(xiàn)象的發(fā)展過程示意圖[76]

上面DBM模型使用的均是單步反應(yīng)模型。2020年，林傳棟和羅開紅[77]利用含化學(xué)反應(yīng)的DBM研究了考慮非平衡效應(yīng)的非穩(wěn)定爆轟現(xiàn)象，該研究中使用了以下兩步鏈?zhǔn)交瘜W(xué)反應(yīng)模型[78]：

(17)

λ′=(1-W)kR(1-λ)exp(-ERT-1)

(18)

式中：ξ′和λ′分別為點(diǎn)火階段以及化學(xué)反應(yīng)放熱階段中的反應(yīng)進(jìn)程參數(shù)隨時(shí)間的導(dǎo)數(shù)；Ts為初始

沖擊溫度；EI為描述化學(xué)誘導(dǎo)過程中溫度敏感性的全局活化能；kI為點(diǎn)火過程方程中指數(shù)前的參數(shù)；W=W(ξ)為階梯函數(shù)，當(dāng)ξ<1時(shí)，W=1,當(dāng)ξ≥1時(shí)，W=0；ER為活化能；kR為放熱過程方程中指數(shù)前的參數(shù)。在該工作中分別研究了擾動(dòng)振幅、波長(zhǎng)以及化學(xué)反應(yīng)放熱對(duì)非穩(wěn)態(tài)爆轟的影響，發(fā)現(xiàn)初始擾動(dòng)幅度僅影響初始階段非穩(wěn)定自持爆轟的形成，當(dāng)初始擾動(dòng)具有較小的波長(zhǎng)時(shí)，壓強(qiáng)在早期階段以更高的振蕩頻率更快地增加，但之后很快減小，并在后期階段變得更小，此時(shí)全局非平衡強(qiáng)度更大，但如果波長(zhǎng)足夠小則全局非平衡強(qiáng)度較小，在這種情況下，最大壓強(qiáng)則展示出相對(duì)小振幅、小平均值以及一個(gè)長(zhǎng)的振蕩周期。此外，還發(fā)現(xiàn)隨著化學(xué)反應(yīng)放熱的增加，壓強(qiáng)和它的振幅增大，非平衡效應(yīng)也增強(qiáng)，但振蕩周期減小。如果擾動(dòng)的波長(zhǎng)足夠小，則不存在橫波或胞格結(jié)構(gòu)，并且二維非定常爆轟會(huì)減弱為一維爆轟。

關(guān)于沖擊波附近的流體動(dòng)力學(xué)和熱力學(xué)非平衡效應(yīng)，林傳棟等[79]通過定義絕對(duì)和相對(duì)偏差度來描述流體系統(tǒng)偏離平衡態(tài)的程度，研究了沖擊波局部和整體的非平衡效應(yīng)以及非組織能量流及其通量，并對(duì)馳豫時(shí)間、馬赫數(shù)、熱導(dǎo)率、黏性和比熱比對(duì)沖擊波處非平衡效應(yīng)的影響進(jìn)行了研究。2021年，吉雨等[80]提出了一個(gè)含化學(xué)反應(yīng)的三維多馳豫DBM，該模型可以自由調(diào)節(jié)普朗特?cái)?shù)和比熱比。通過引入Arrhenius不可逆、單步化學(xué)反應(yīng)模型，模擬自由下落過程中的化學(xué)反應(yīng)、Couette流、一維穩(wěn)態(tài)和非穩(wěn)態(tài)爆轟以及在封閉立方體中的三維球形爆炸驗(yàn)證了新模型的正確性。圖10給出該模型的一組模擬結(jié)果：外爆過程中3個(gè)時(shí)刻的壓強(qiáng)場(chǎng)分布。

圖10 外爆過程中壓強(qiáng)場(chǎng)的演化過程[80]

4 結(jié)論與說明

超聲速、高瞬變等極端條件引發(fā)一系列值得深入研究的可壓、非平衡復(fù)雜流動(dòng)行為；非線性效應(yīng)、離散效應(yīng)、強(qiáng)耦合效應(yīng)是其重要特征。隨著流體行為非平衡程度增強(qiáng)，只關(guān)注分布函數(shù)守恒矩(密度、動(dòng)量和能量)演化的描述方法在物理功能方面越來越不能滿足需求。其表現(xiàn)之一是：不考慮非平衡效應(yīng)或者非平衡效應(yīng)處理不當(dāng)，直接影響著密度、流速、溫度、壓強(qiáng)這些常用宏觀量結(jié)果的精度。另外，這些以前因模型、方法原因而知之甚少的非平衡行為特征蘊(yùn)含著大量待開發(fā)的物理功能。發(fā)展DBM等介尺度建模和模擬方法，研究這些非平衡行為特征；基于獲得的新認(rèn)識(shí)，開發(fā)新的物理功能正在成為該研究方向的重要內(nèi)容。

需要說明的是：① DBM建模方法主要是針對(duì)宏觀連續(xù)建模失效或物理功能不足、而微觀分子動(dòng)力學(xué)模擬又因?yàn)檫m用尺度問題無能為力的“介尺度”情形，應(yīng)“介尺度”非平衡系統(tǒng)研究需求而設(shè)計(jì)的物理模型構(gòu)建方法；② DBM提供的物理信息量介于宏觀連續(xù)描述和微觀分子動(dòng)力學(xué)之間；③ 相對(duì)于宏觀描述，DBM從一個(gè)更寬的視角來觀測(cè)系統(tǒng)；④ DBM中非守恒矩描述的必要性和收益均隨著非平衡程度增高而增加；⑤ 與不含化學(xué)反應(yīng)的復(fù)雜流動(dòng)相比，燃燒系統(tǒng)的熱力學(xué)非平衡多了一個(gè)來源：快速反應(yīng)。當(dāng)化學(xué)反應(yīng)速率加快到一定程度，便不能再假設(shè)在化學(xué)反應(yīng)進(jìn)程的每一步系統(tǒng)都回到了熱力學(xué)平衡態(tài)；⑥ 目前已公開發(fā)表的DBM研究涉及的只是相對(duì)稀薄和快速流動(dòng)引發(fā)的熱力學(xué)非平衡。同時(shí)，盡管DBM可以在非平衡行為描述的廣度與深度2個(gè)方面超越Navier-Stokes方程，但鑒于研究的階段性，目前已公開發(fā)表的燃燒系統(tǒng)的DBM研究主要集中在DBM的描述廣度方面；且其應(yīng)用研究也遠(yuǎn)不如在相分離[28,81-84]和流體不穩(wěn)定性研究[21-22,24,26,64,85-90]方面充分，除此之外，目前的DBM燃燒多相模擬中的多相主要指燃燒過程中的多組分流體。因而，包含化學(xué)反應(yīng)非平衡的DBM建模與模擬、在非平衡深度描述方面超越Navier-Stokes方程的DBM建模、分布函數(shù)f對(duì)平衡態(tài)feq的偏離(f-feq)對(duì)外力項(xiàng)效應(yīng)影響較大情形的DBM建模與模擬、多個(gè)相(態(tài))的DBM燃燒模擬、燃燒多相流系統(tǒng)的深度非平衡動(dòng)理學(xué)機(jī)理與應(yīng)用均是下一步研究工作的重點(diǎn)。

致 謝

感謝閆鉑、張玉東、賴惠林、羅開紅、陳正、王健平、王兵、聶百勝、王成、孫遠(yuǎn)翔等在燃燒建模與模擬方面和王立鋒、趙英奎等在內(nèi)爆動(dòng)理學(xué)等方面的有益討論。