低頻爆轟不穩定性形成機理的數值模擬研究*

張允禎，程 杪，榮光耀，王健平

（北京大學工學院力學與工程科學系，北京 100871）

爆轟是一種激波與波后化學反應區強烈耦合的物理現象。爆轟現象自被發現以來，已受到了廣泛的關注，并已形成了相關的研究理論。目前而言，被普遍認可并被廣泛應用的爆轟理論主要有兩種，分別為C-J理論和ZND模型[1-2]。由于爆轟相較于燃燒具有熵增小、熱效率高等特點[3]，因此研究人員也一直致力于將爆轟引入到化學推進系統中，以期實現發動機性能的大幅提升。而近年來，連續爆轟發動機（rotating detonation engine）更是成為國際上研究的熱點。

圖1為連續爆轟發動機的原理示意圖。在發動機運行過程中，在頭部進氣壁面附近有一道或多道沿圓周方向連續旋轉傳播的爆轟波頭（detonation front），推進劑從頭部進氣壁面連續不斷地噴入燃燒室內，在爆轟波頭前形成一個新鮮氣體層（fresh gas layer）。在爆轟波頭的斜后方附著有一道斜激波（oblique shock wave）和一道接觸間斷（contact surface），接觸間斷兩側分別為上周和本周的爆轟產物（burnt gas），高溫高壓的爆轟產物從尾部高速噴出，從而獲得所需的推力。

圖1 連續爆轟發動機原理（二維圓柱流場及其展開結構）Fig.1 Schemetic of RDE(2D cylindrical flow field and its unfolding structure)

20世紀50年代末至60年代，Voitsekhovskii[4]和Nicholls等[5]分別對連續爆轟發動機進行了驗證性研究。此后由于研究條件的限制，連續爆轟發動機的相關研究陷入近乎停滯的狀態，直到Bykovskii等[6]重新在實驗上實現了較長時間的連續旋轉爆轟。此后世界各地學者們[7-17]對相關領域開展了大量的理論和實驗研究。

對于一種推進系統，其能夠穩定可靠地運行是研究人員追求的最終目標，因此連續爆轟發動機的穩定性研究是一個不可避免的課題。在以往針對連續爆轟發動機的研究中，人們發現，在連續爆轟發動機運行過程中存在各式各樣的不穩定性，并針對這些不穩定性開展了一定的研究。在這個過程中，一些研究人員對這些不穩定性現象進行了分類。Wang等[18]根據實驗測得的不穩定信號的振蕩頻率將這些不穩定性分為低頻、中頻、高頻不穩定性，并指出低頻爆轟不穩定性是由燃燒室和進氣室的聲學耦合引起的。Anand等[19-21]按照壓力信號的特征將所發現的不穩定現象分為4類，并針對這些不穩定性開展了持續的、系統的實驗研究。針對低頻爆轟不穩定性（low frequency instability,LFI），Anand等[21]通過實驗開展了深入研究，他們認為產生這種低頻爆轟不穩定性的一個原因是由于燃燒室頭部的激波與爆轟波相互作用所致。此外，Zhang等[22]的研究表明，燃燒室頭部的激波會引起新鮮氣體層的不規則分布，進而導致了低頻爆轟不穩定性，但是研究未能揭示不規則新鮮氣體層與壓力峰值低頻振蕩之間的直接關聯。此外，實驗和數值上均未給出激波的來源、發展及導致低頻爆轟不穩定性的具體作用過程。本文將通過數值模擬的方式針對低頻爆轟不穩定性現象開展進一步的研究。

1 數值方法和物理模型

1.1 數值方法

本文采用H2/Air 9組分19步基元化學反應模型針對同軸圓環腔燃燒室開展二維的數值模擬研究。爆轟波的傳播速度非常快，一般都在千米每秒的量級，因此在對爆轟波的數值模擬中，一般忽略組分的擴散、黏性以及熱傳導[23]的影響，采用二維圓柱坐標系下的帶源項的守恒型的Euler方程作為控制方程[22,24]：

其中

式中：U為守恒量，F和G為通量向量，S為化學反應源項；ρ為混合氣體總密度，ρi為各組分氣體的密度，需滿足歸一化條件；E為混合氣體單位體積的總能，h為混合氣體總比焓；uθ為周向速度分量，uz為軸向速度分量；為第i個組分的化學生成速率，由基元化學反應模型得到；hi為第i個組分的比焓，與混合氣體的溫度有關。由于溫度、壓力、能量強烈耦合，因此聯立式(3)、(5)，采用Newton-Rapson迭代法求解溫度。熱力學參數采用溫度的多項式進行擬合。

在空間上采用Steger-Warming矢通量分裂法對通量向量進行分裂，采用五階WENO格式對通量進行重構，時間上采用2階TVD Runge-Kutta方法推進求解。在對化學反應的處理上，采用解耦算法將流動和化學反應進行解耦，并采用劉君等[25]提出的半隱式格式求解化學反應源項以解決求解源項時產生的剛性問題，并同時提高計算效率。

1.2 物理模型

本文將對同軸圓環腔燃燒室內存在的低頻爆轟不穩定性開展數值模擬研究。同軸圓環腔燃燒室的環腔內外徑之差遠小于半徑，因此在計算中可以將其簡化為一個沒有厚度的圓柱面，并忽略徑向的影響，于是可將此圓柱面沿母線展開，在二維平面上進行計算。圖2為一軸向長度（L）為0.04 m、半徑（R）為0.01 m的燃燒室的平面展開計算域。其中下邊界為進氣壁面，上邊界為燃燒室尾部出口。按化學恰當比預混好的H2/Air混合氣從進氣壁面噴入燃燒室，入流總壓ptot=5×105Pa，入流總溫Ttot=360 K。

圖2 初始時刻計算域Fig.2 Computational domain at the initial time

初始時刻，燃燒室內充滿按照化學恰當比預混好的新鮮可燃氣體，初始壓強p0和初始溫度T0分別為1×105Pa和300 K。在頭部壁面附近覆蓋一個典型的一維爆轟波，用以起爆周向傳播的旋轉爆轟波。

邊界條件采用如下設置。

（1）下邊界為入流邊界，采用收斂噴管入流邊界。將每個網格點視為一個收斂噴管，在每個網格點上對物理量按照收斂噴管的條件進行設置。進氣壁面上的壓強pw為收斂噴管的背壓，入流總壓為ptot，入流總溫為Ttot。收斂噴管的計算中有兩個關鍵參數，入流總壓ptot和臨界壓強pcr：(a)若pw＞ptot，新鮮氣體無法噴入燃燒室，此時采用固壁反射邊界條件，速度矢量在壁面兩側沿法向相反，沿切向相等，密度、壓力、溫度等標量值在壁面兩側對稱相等；(b)若pw＞ptot＞pcr，新鮮氣體以亞聲速噴入燃燒室；(c)若pw＜pcr，噴管壅塞，新鮮氣體以聲速噴入燃燒室。

（2）上邊界為尾部出流邊界，采用亞聲速/超聲速出流邊界條件。(a)若出口處的馬赫數Ma＜1，氣體以亞音速出流，邊界處的壓力取為周圍環境的壓力，其余物理量由邊界內部的網格點上對應的物理量插值得到；(b)若Ma＞1，氣體以超聲速噴出燃燒室，下游的流場不能影響上游的流場，因此虛擬網格的值可隨意設定，一般設定與出口處的值相等。

（3）左右邊界為圓柱壁面上同一條母線，設定周期性邊界條件，即將一側的邊界內側節點上的值賦值給另一側的邊界外側的虛擬網格點上。

由于采用基元化學反應模型，計算量較大，因此采用MPI并行計算，以提高計算效率。圖3為采用0.1 mm和0.2 mm網格計算所得流場的溫度云√圖和特征壓力梯度對數云圖。特征壓力梯度對數定義為[26]：dp?=0.8Exp[?100|?p|/max|?p|]，其中|?p|=(?p/?x)2+(?p/?y)2。從溫度云圖可以看到，兩種尺寸的網格均能夠模擬出在燃燒室頭部周向傳播的旋轉爆轟波，且二者得到的流場基本一致。但從特征壓力梯度對數云圖來看，0.2 mm的網格得到的流場的波系結構較粗糙，不如0.1 mm的網格得到的流場展示出來的細節豐富。由于本文的研究中需要針對流場中的細微結構進行分析，因此采用0.1 mm的網格進行計算。

圖3 不同網格尺寸的計算結果（溫度云圖和特征壓力梯度對數云圖）Fig.3 Calculation results of different grid sizes(temperature nephogram and pressure gradient logarithmic nephogram)

2 結果與討論

在上述的初始及邊界條件下，最終在燃燒室內形成了一個沿周向傳播的爆轟波，記錄燃燒室進氣壁面上某采樣點的壓強隨時間的變化，得到p-t曲線，如圖4(a)所示。圖4(b)為Anand等[21]的H2/Air連續爆轟發動機實驗結果，3種不同顏色為沿環向均布的3個位置處的壓力時程曲線。數值模擬結果與實驗結果高度類似，爆轟波峰值均呈現出周期性起伏振蕩的現象。計算爆轟波在燃燒室內傳播的頻率fD（會受到燃燒室尺寸的影響）和壓力振蕩頻率fosc，數值模擬的結果分別為fD≈30 kHz和fosc≈1.5 kHz，兩者之比fosc/fD≈0.05；Anand的實驗結果兩者分別約為2 kHz和150 Hz，兩者之比。可以看到，數值和實驗結果均非常小且十分相近，在同一量級。因此，一般將這種現象稱為低頻爆轟不穩定性。

圖4 采樣點（傳感器處）的壓強-時間曲線Fig.4 The pressure-time track of sampling points(or the points where the sensors are placed)

2.1 爆轟波頭上壓強的周期性變化

為探究這種低頻爆轟不穩定性產生的原因，針對燃燒室內爆轟波的傳播過程進行重點研究。圖5為340～402 μs進氣壁面上最大壓強值隨時間的變化。可以發現，進氣壁面上的最大壓強也呈現出了周期性起伏振蕩的現象。選取黑色虛線框內（374～385 μs）的一個周期，詳細觀察燃燒室內的流場。圖6為這段時間內爆轟波頭在燃燒室頭部的傳播過程，其中白色曲線為接觸間斷，即已燃氣與未燃氣的分界線（因此白色曲線與進氣壁面之間即為新鮮氣體層）。圖7為這段時間內爆轟波頭上的壓強沿軸向的分布情況，橫軸代表軸向的距離，縱軸代表相應距離上爆轟波頭上的壓強。在這種工況下，一般認為壓強超過4.5 MPa為爆轟波。

圖5 進氣壁面上最大壓強隨時間的變化Fig.5 Track of the peak pressure on the inlet wall

圖6 374～384 μs爆轟波頭的傳播情況Fig.6 Propagation of detonation front between 374?384 μs

從圖6和圖7可以看到，爆轟波頭上的壓強分布在時間上并不一致，會存在一些向進氣壁面運動的高壓區。在374 μs時，爆轟波頭的高壓區正與進氣壁面接觸（圖6），因此此時進氣壁面最大壓強達到最大值（圖5）。在374～378 μs之間，爆轟波越過新鮮氣體層上的凹點（以下將之稱為進氣阻滯點(injet blocking point,IBP)，在2.3節中會解釋原因），波前的新鮮氣體層的寬度急劇變寬（圖6）。此時，在爆轟波頭與新鮮氣/已燃氣分界線相交的地方產生了一道局部爆炸。這個現象可以直觀地從圖7中377～379 μs的波頭壓強分布曲線上看出。隨后，這道局部爆炸產生的環形激波使得爆轟波頭上產生一個高壓區，這個高壓區向進氣壁面運動，與進氣壁面相遇（如圖7所示）從而使得進氣壁面上的最大壓強值起伏變化。由于爆轟波周期性與進氣阻滯點相遇，爆轟波頭上壓強的軸向分布周期性變化，從而進氣壁面的最大壓強便產生了周期性的振蕩。

圖7 374～385 μs爆轟波頭上的壓強軸向分布Fig.7 Axial distribution of the pressure on the detonation front between 374?385 μs

像這種爆轟波頭上產生局部爆炸的現象，其他學者在研究中也有發現。Hishida等[27]在他們的數值研究中曾經發現過類似的現象并采用了較精細的網格進行了研究，認為爆轟波越過新鮮氣體層上的凹點時，爆轟波在新鮮氣/已燃氣分界線處會產生未燃氣包，從而會形成局部爆炸。最近，Athmanathan等[28]在開展的連續爆轟發動機可視化實驗中觀察到，爆轟波頭上存在所謂的強爆轟區，且這些強爆轟區會向燃燒室的頭部方向運動。他們認為這種現象是由于爆轟波頭前新鮮氣體層高度的變化導致的。

2.2 低頻爆轟不穩定性產生的機理

圖8為340～402 μs之間進氣壁面上最大壓強（爆轟波頭與進氣壁面相交處的壓強）的變化，藍色虛線為爆轟波頭經過進氣阻滯點的時刻。將爆轟波經過兩個進氣阻滯點之間的時間記為 ?t，相鄰的兩個進氣阻滯點間的距離為 ?x，爆轟波速D視為恒定的。由以上研究可以看到，進氣壁面上相鄰兩個進氣阻滯點間的距離 ?x基本一致且保持不變，因此爆轟波經過兩個進氣阻滯點間的時間 ?t=?x/D也基本恒定，這一點可以很直觀地從圖8中看出。可以看到，爆轟波每經過任意相鄰的兩個進氣阻滯點的時間里，進氣壁面上最大壓強的歷程基本一致，總體上呈現出較規則的正弦振蕩。

圖8 進氣壁面上最大壓強的變化（數值模擬結果）Fig.8 Track of the peak pressure on the inlet wall(numerical result)

如圖9(a)所示，當爆轟波處于兩進氣阻滯點（紅色方塊處）之間的某個位置時，設爆轟波頭距離上個進氣阻滯點（設為第n個）的距離為l。此時，爆轟波距離經過第n個進氣阻滯點時已經過了t′=l/D的時間，于是此時進氣壁面上的最大壓強即為如圖9(b)所示的 (n?t+t′)時刻黑色虛線所指示的壓強。

圖9 采樣點壓強振蕩的形成機理Fig.9 The schematic diagram of mechanism of the pressure oscillation at the sampling point

若此時爆轟波頭正好與進氣壁面上的采樣點（sampling point，SP，圖9(a)中為黑色方塊處）相遇，則此刻采樣點的壓強正好達到壓強峰值，且該壓強峰值的大小為圖9(b)所示的n?t+t′時刻黑色虛線所指示的壓強。在進一步的研究中發現，這些進氣阻滯點產生的位置會在進氣壁面上沿圓周方向緩慢移動，如圖10所示，其中w(H2)為H2的質量分數。因此對于進氣壁面上某固定的采樣點來說，爆轟波每次與之相遇時，爆轟波頭相對于上一個進氣阻滯點的位置l會緩慢變化。對于每個確定的l(0≤l≤?x)，都能夠在圖9(b)上有一個對應的壓強，此壓強即為本周爆轟波旋轉傳播過程中采樣點的壓強峰值。隨著l從 ?x逐漸減小到零，采樣點的壓強峰值便完成了一個周期的振蕩。又由于進氣阻滯點所產生的位置的移動速度v非常小，因此采樣點的壓強-時間曲線上壓強峰值的振蕩頻率f=v/?x也就很低，從而產生了所謂的低頻爆轟不穩定性。

圖10 爆轟波每旋轉兩周進氣阻滯點產生的位置Fig.10 The position of IBP in every two cycles of detonation wave rotation

2.3 不規則新鮮氣體層的形成

從上述研究中可以看到，新鮮氣體層的不規則分布對于低頻爆轟不穩定性的產生起著關鍵作用。本節將闡釋不規則新鮮氣體層的形成原因。

圖11為進氣壁面上進氣阻滯點產生時該點周圍H2質量分數的分布云圖。從圖11中可以看出，391～393 μs之間，進氣壁面上Y=0.03 m附近的進氣被阻滯，該區域兩側的進氣均先于該區域的進氣，其中0.029 m處的進氣被阻滯的時間最長。分別在Y=0.027,0.029,0.032 m處設置采樣點，記錄進氣被阻滯處前后的壓強變化，并將之與穩定的進氣情況進行對比，如圖12所示。進氣壁面壓強低于進氣總壓（紅色虛線，500 kPa）被認為氣體可以噴入燃燒室。

圖11 進氣阻滯點的產生Fig.11 Generation of injet blocking points

圖12 進氣被阻滯處及左右各一點在進氣阻滯點產生時的壓強變化和穩定進氣情況Fig.12 Pressure track of the place where intake process is interrupted and one point on each side when IBP is being generated,and the stable intake situation

與穩定情況進行對比，可以發現，發生進氣阻滯的情況下，開始進氣的前后進氣壁面上的壓強在時間上分布并不平滑，有一系列的激波沿著與爆轟波相反的方向從0.032 m處向0.029 m處運動（圖12中綠色圓圈所示）。這些激波傳到0.032 m處時，該點的壓強還遠高于進氣總壓，因此只是在壓強曲線位于進氣總壓上面的部分引起了一定的振蕩，在一定程度上使得該區域的進氣過程滯后了一些，但對進氣的阻滯效果十分有限；當這些激波傳播到0.029 m處時，該點的壓強剛好下降到進氣總壓，新鮮氣體即將噴入燃燒室，但是這些激波的到來引起了即將進氣的地方壓強產生強烈振蕩，并使壓強保持在進氣總壓以上很長一段時間，即將噴入的新鮮氣體被阻止進入燃燒室，進氣被強行打斷，于是在391～393 μs形成了圖11中0.029 m處的進氣阻滯點；同時，當這些激波經過0.029 m并引起進氣阻滯的時候，上游（0.032 m附近）的壓強卻已經下降到進氣總壓以下，新鮮氣體開始噴入燃燒室；隨后，當這些激波越過0.029 m（進氣被長時間阻滯處）到達下游的0.027 m處時，如前所述，該點的壓強也早已下降到了進氣總壓以下并開始進氣，并且這些激波隨著傳播距離的增加其強度也逐漸衰減，因此在0.027 m處不能阻礙進氣過程，只是在該點的壓強曲線上位于進氣總壓以下的部分引起了一定的壓強回升。因此，這些激波僅在與新鮮氣體層頂點（即壓強剛下降到進氣總壓以下即將開始進氣的點）相遇的地方間斷性地導致了進氣阻滯點的產生。

接下來，將從流場內波系結構的發展出發，揭示這些激波的來源、發展與導致進氣阻滯點產生的詳細過程。首先，由Chen等[29]的研究可知，當反傳激波與爆轟波相遇后，反傳激波的強度將得到大幅增強。圖13為362～393 μs燃燒室頭部波系結構的發展過程。從圖13(a)中可以看到，爆轟波后存在一些強度較強的激波（黑色箭頭所示），這些激波不斷地向后傳播，其強度不斷衰減，最終進入新鮮氣體層（371～374 μs），在下周的爆轟波前形成了一個橫波區。在381～384 μs之間（圖13(b)），這些波前橫波（黑色箭頭所示）與爆轟波相遇后被大幅增強，在波后重新形成了一個更強的橫波區（紅色箭頭所示）。之后，在385～393 μs之間（圖13(c)），這些被增強后的橫波（黑色箭頭所示）向進氣壁面運動，與進氣壁面相撞，發生反射。反射點與新鮮氣體層的頂點（白色箭頭所示）相遇，使得該點處的壓強長時間維持在進氣總壓以上，該點的進氣過程被長時間打斷，于是形成了新鮮氣體層上的進氣阻滯點。圖13(c)的白色圓圈內可以看到一道較強的激波在進氣壁面上的反射點與新鮮氣體層頂點相遇的情景。

圖13 燃燒室頭部壓強對數云圖Fig.13 Logarithmic nephogram of pressure at the head of combustor

這些激波在燃燒室頭部不斷地與爆轟波相向傳播，其強度在傳播過程中衰減，而在與爆轟波相遇之后又被增強。這個過程循環往復，于是這些激波得以持續存在于燃燒室中，使進氣壁面上不斷產生進氣阻滯點，從而導致了新鮮氣體層的不規則分布，進而導致了低頻爆轟不穩定性的產生。

3 結 論

本文通過二維數值模擬對連續爆轟發動機中常見的低頻爆轟不穩定性進行了研究，揭示了低頻爆轟不穩定性形成的機制及詳細過程。

（1）燃燒室頭部附近存在一些與爆轟波傳播方向相反的反傳激波，這些激波與爆轟波相遇后會被增強，從而能夠持續存在于燃燒室中。這些激波會與進氣壁面發生反射，反射點與新鮮氣體層的頂點相遇后會阻斷該點的進氣過程，從而在該點產生進氣阻滯點，導致新鮮氣體層的不規則分布。

（2）新鮮氣體層的不規則分布會導致爆轟波頭上周期性形成高壓區，導致爆轟波頭上的壓強隨著進氣阻滯點的分布而產生周期性變化，進而引起進氣壁面上最大壓強的起伏振蕩。

（3）隨著這些進氣阻滯點的生成位置沿圓周方向緩慢地移動，爆轟波每旋轉一周經過采樣點時，距離上一個進氣阻滯點的距離l會緩慢變化，波頭上與進氣壁面相接觸的位置的壓強（即采樣點的壓強，亦即距離l所對應的壓強）便產生了低頻率的周期性起伏振蕩，即形成了低頻爆轟不穩定性。

感謝國家超級計算天津中心的大力支持。