燃燒多相流的介尺度動(dòng)理學(xué)建模研究進(jìn)展

許愛國，單奕銘, 陳鋒，甘延標(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)步，人類飛行器的速度越來越高。“高速氣流中的點(diǎn)火與持續(xù)燃燒成為空天推進(jìn)中長久的關(guān)鍵問題”[1]。這些關(guān)鍵問題，往往涉及極端條件。極端條件，是相對于常溫常壓常速等普通條件而言的，它至少在某一方面嚴(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ì)對飛行器表面產(chǎn)生燒蝕作用，使其防熱結(jié)構(gòu)遭到破壞，從而改變飛行器氣動(dòng)性能，影響其穩(wěn)定性和飛行安全性。并且在高溫下飛行器表面的防熱材料也會(huì)影響邊界層內(nèi)氣體化學(xué)反應(yīng)的催化效應(yīng)，不同壁面催化條件會(huì)對氣動(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)，這就需要對飛行進(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)，飛行器的外流場表現(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ù)雜物理場分析角度，介紹離散玻爾茲曼建模方法(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 宏觀與微觀建模

人們先前對燃燒和爆轟問題的研究主要依賴于實(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]，其中分子間作用勢的確立是模型構(gòu)建的關(guān)鍵，分子間作用勢有效半徑的選取是模擬過程中的關(guān)鍵，有效半徑的選取通常以滿足需求的最小值為最佳原則。微觀研究可以幫助建立反應(yīng)速率方程，提供完備信息，但由于運(yùn)算量太大，對內(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方程對應(yīng)的宏觀流體方程組就是Euler方程組；如果系統(tǒng)時(shí)刻處于熱動(dòng)平衡態(tài)附近，二階及以上非平衡效應(yīng)很弱以至于可以忽略，只需考慮一階非平衡效應(yīng)，那么此時(shí) Boltzmann方程對應(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)”(相對于所關(guān)注的結(jié)構(gòu)尺度，平均分子間距不再是可以忽略的小量，不能再使用連續(xù)性假設(shè))，模式快到一定程度會(huì)引發(fā)“熱力學(xué)非平衡效應(yīng)”(相對于所關(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)槟Ｐ秃头椒ǖ娜狈?不成熟)而研究薄弱。這些薄弱不僅在一定程度上阻斷了對微觀到宏觀之間關(guān)聯(lián)的認(rèn)識(shí)，也在一定程度上阻礙了本構(gòu)模型的機(jī)理化、科學(xué)化；同時(shí)，微流控等技術(shù)的快速發(fā)展在提醒人們，這些(相對于宏觀)特征更加豐富但以前知之甚少的“介尺度”非平衡、不連續(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é)矩在模型簡化過程中保值；以該組動(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)心的物理特征量使用簡化前和簡化后的模型計(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,…]對應(yīng)研究中關(guān)心的、建模過程中要保值的動(dòng)理學(xué)矩，v為分子速度；fi=f(vi)，vi為第i個(gè)離散速度。第3步，是DBM建模的目的和核心，給出非平衡狀態(tài)和行為描述的具體方案。DBM建模主要針對的是宏觀連續(xù)建模失效或物理功能不足、而微觀分子動(dòng)力學(xué)模擬又因?yàn)檫m用尺度問題無能為力的“介尺度”“兩難”情形。DBM提供的物理信息量介于宏觀連續(xù)描述和微觀分子動(dòng)力學(xué)之間。相對于宏觀描述，DBM從一個(gè)更寬的視角來觀測系統(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)常被誤解為LBM方程解法，LBM建模方法逐漸改稱為LB動(dòng)理學(xué)建模(LB Kinetic Modeling, LBKM)、DBM。圖1中DBM建模與分析方法軸上2012以后的點(diǎn)對應(yīng)燃燒系統(tǒng)DBM建模文章發(fā)表的年份。需要指出的是，統(tǒng)計(jì)物理學(xué)領(lǐng)域作為粗粒化建模的格子氣方法，不會(huì)因?yàn)镈BM或其他建模方法的出現(xiàn)而消失，它的思路永遠(yuǎn)相對獨(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é)勢引起的。在已發(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)勢也就越明顯。

3 部分研究結(jié)果

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

2014年，林傳棟等[70]基于極坐標(biāo)系構(gòu)建了一個(gè)用于模擬燃燒的DBM，模擬了典型的內(nèi)爆和外爆過程，觀測并分析了化學(xué)反應(yīng)熱、物質(zhì)輸運(yùn)、熱傳導(dǎo)與幾何(匯聚或發(fā)散)效應(yīng)的合作與競爭。在外爆過程中觀測到了與一維爆轟情形不同的、幾何效應(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)勢, 則會(huì)出現(xiàn)熄爆現(xiàn)象。隨著爆轟波向外傳播, 幾何發(fā)散效應(yīng)逐漸減弱, 發(fā)生熄爆的可能性降低。發(fā)現(xiàn)在高度對稱的系統(tǒng)中，球心處系統(tǒng)始終處于熱力學(xué)平衡態(tài)。2018年，許愛國等提出一個(gè)單弛豫時(shí)間球坐標(biāo)DBM模型[71]。在內(nèi)爆與外爆過程中，幾何匯聚與發(fā)散效應(yīng)起到一個(gè)“外場力”的作用，更加細(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í)觀測反應(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ò)散方程等)均是該模型的特例，借助該模型可以很方便地觀測和研究與宏觀行為相伴隨的各種熱力學(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)行比對，定壓燃燒的火焰溫度(T)與實(shí)驗(yàn)結(jié)果符合較好[73]，如圖3所示。因?yàn)椴煌姆磻?yīng)物得以分別描述，所以該模型既可以描述預(yù)混燃燒，又可以描述非預(yù)混燃燒。

2.7 準(zhǔn)確度考察（加樣回收率試驗(yàn)） 取錐形瓶6個(gè)，分別精密加入新配制的大黃素-8-O-β-D-葡萄糖苷對照品溶液（66 μg∕mL）6.36 mL、大黃素甲醚-8-O-β-D-葡萄糖苷對照品溶液（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)樣測定，按下式計(jì)算回收率。

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

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

為了估算爆轟波波峰的相對高度，定義了峰高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ù)對爆轟演化過程的影響。使用的化學(xué)反應(yīng)率模型為

(12)

其中:

(13)

(14)

(15)

式中：Tth為起爆溫度；k為化學(xué)反應(yīng)速率系數(shù)；h1和h2分別代表了k的峰值和谷值；T1和T2分別是h1和h2對應(yīng)的溫度。研究了圖7(a)～圖7(d)所示4種反應(yīng)速率系數(shù)下的爆轟現(xiàn)象。通過比較可以發(fā)現(xiàn)，對于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)量流以及無組織能量流。對于系統(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ù)對動(dòng)力學(xué)量的作用是降低馮·紐曼峰的高度(減小密度、壓強(qiáng)和速度)，加寬反應(yīng)區(qū)，抑制化學(xué)反應(yīng)進(jìn)程。

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

2019年，張玉東等[76]提出了一個(gè)用于模擬爆轟的一維離散玻爾茲曼模型。通過對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)象。對于反常爆轟來說，其并沒有一個(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)，對于常規(guī)的爆轟現(xiàn)象，化學(xué)反應(yīng)主要發(fā)生在S1和S2階段中，但對于該反常爆轟，在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)振幅、波長以及化學(xué)反應(yīng)放熱對非穩(wěn)態(tài)爆轟的影響，發(fā)現(xiàn)初始擾動(dòng)幅度僅影響初始階段非穩(wěn)定自持爆轟的形成，當(dāng)初始擾動(dòng)具有較小的波長時(shí)，壓強(qiáng)在早期階段以更高的振蕩頻率更快地增加，但之后很快減小，并在后期階段變得更小，此時(shí)全局非平衡強(qiáng)度更大，但如果波長足夠小則全局非平衡強(qiáng)度較小，在這種情況下，最大壓強(qiáng)則展示出相對小振幅、小平均值以及一個(gè)長的振蕩周期。此外，還發(fā)現(xiàn)隨著化學(xué)反應(yīng)放熱的增加，壓強(qiáng)和它的振幅增大，非平衡效應(yīng)也增強(qiáng)，但振蕩周期減小。如果擾動(dòng)的波長足夠小，則不存在橫波或胞格結(jié)構(gòu)，并且二維非定常爆轟會(huì)減弱為一維爆轟。

關(guān)于沖擊波附近的流體動(dòng)力學(xué)和熱力學(xué)非平衡效應(yīng)，林傳棟等[79]通過定義絕對和相對偏差度來描述流體系統(tǒng)偏離平衡態(tài)的程度，研究了沖擊波局部和整體的非平衡效應(yīng)以及非組織能量流及其通量，并對馳豫時(shí)間、馬赫數(shù)、熱導(dǎo)率、黏性和比熱比對沖擊波處非平衡效應(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)場分布。

圖10 外爆過程中壓強(qiá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建模方法主要是針對宏觀連續(xù)建模失效或物理功能不足、而微觀分子動(dòng)力學(xué)模擬又因?yàn)檫m用尺度問題無能為力的“介尺度”情形，應(yīng)“介尺度”非平衡系統(tǒng)研究需求而設(shè)計(jì)的物理模型構(gòu)建方法；② DBM提供的物理信息量介于宏觀連續(xù)描述和微觀分子動(dòng)力學(xué)之間；③ 相對于宏觀描述，DBM從一個(gè)更寬的視角來觀測系統(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研究涉及的只是相對稀薄和快速流動(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對平衡態(tài)feq的偏離(f-feq)對外力項(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é)等方面的有益討論。