同軸離心式噴嘴熱聲不穩定性遞歸分析

楊向明，楊尚榮，楊岸龍，艾春安

(1. 火箭軍工程大學， 西安 710025； 2. 西安航天動力研究所液體火箭發動機技術重點實驗室， 西安 710100)

0 引 言

熱聲不穩定是火焰、流場和聲場間非線性耦合的結果，其本質是非線性系統[1]。因此，可以采用非線性時間序列方法，從動力系統角度分析非線性熱聲系統的動力學特征，尋求不同動力學狀態間轉變的控制參數，進而對熱聲不穩定進行提前預報和控制[2]。

非線性時間序列分析方法中，非線性不變量，如分形維數、最大Lyapunov指數分別表征系統在相空間中的靜態和動態結構，熵代表信息生成率，均被用來分析熱聲系統的動力學特性。Kabiraj等[3]利用關聯維數和最大Lyapunov指數研究了管道層流預混燃燒系統的熱聲振蕩現象。Gotoda等研究了貧燃預混旋流穩定湍流火焰的動力學特性，發現平移誤差[4]、排列熵[5]均可作為監測貧燃熄滅極限的指標。

另一類基于動力系統遞歸性質的非線性時間序列方法，如遞歸圖和遞歸量化分析方法，可以從熱聲系統中的標量時間序列中獲得更多系統動力學方面的特性，在熱聲系統動力學分析中得到了較多應用。文獻[6-7]采用遞歸圖分析方法，在鈍體穩焰湍流燃燒器中，發現在熱聲不穩定發生之前，高幅值的周期振蕩總會陣發性出現，直到熱聲不穩定發生為止，這種陣發性現象也被多項實驗證實[8-9]。Sen等[10]利用遞歸量化方法對管道逆擴散火焰進行了研究，發現了極限環、陣發和同宿異軌等現象。Gotoda等[11]利用遞歸圖、多尺度熵等方法研究了貧燃預混旋流穩定湍流火焰的動力學特性，發現在接近貧燃熄滅邊界時，存在混沌向周期振蕩轉變的隨機過程，繼續增加當量比，周期振蕩通過準周期振蕩又轉變為混沌[12]。Unni等[13]分別利用廣義的赫斯特指數和遞歸量化分析方法，在同一框架內對湍流燃燒系統中的熱聲不穩定和貧燃熄滅[14]兩種現象進行了預報。結合遞歸圖、遞歸網絡方法和同步理論[15-16]，學者發現非穩態熱釋放率對聲壓的影響更大(相比聲壓對非穩態熱釋放率的影響)，因此認為消除熱聲不穩定時控制火焰比控制聲場更有效。

液體火箭發動機中的熱聲不穩定是發動機研制過程中必須解決的關鍵難題。同軸離心式噴嘴作為最常用的液體火箭發動機雙組元噴嘴，對其燃燒穩定性的研究一直受到學者們關注。王楓等[17]利用相似準則設計了縮比燃燒室，發現噴嘴縮進長度對燃燒穩定性裕量有很大影響。王迪等[18]開展了氣液同軸離心式單噴嘴模型發動機試驗，發現噴嘴縮進長度的增加對縱向高頻燃燒不穩定性產生阻尼作用。燃燒室長度對縱向高頻燃燒不穩定性的影響比噴嘴縮進長度更明顯。王延濤等[19]試驗發現同軸離心噴嘴的自發激勵高頻燃燒不穩定過程出現了“滯后”現象，不穩定性工況下的火焰均為脫口火焰，火焰特征長度約等于噴嘴出口到脫口火焰團上沿的距離。

本文開展同軸離心式噴嘴穩定性試驗研究，采用遞歸圖和遞歸量化分析方法對實驗結果進行討論，從動力系統角度分析同軸離心式噴嘴火焰動力學特性，為液體火箭發動機熱聲不穩定性預報和控制提供理論支持。

1 實驗系統和數據處理方法

1.1 實驗裝置

燃燒試驗系統原理如圖1所示。模擬燃燒室為敞口圓筒形結構，直徑Dc= 156 mm，高Hc=277 mm，垂直安放在模擬噴注器面板(水冷平板)上。試驗用噴嘴(見圖2)為帶縮進室的同軸離心噴嘴，由外噴嘴和內噴嘴組成，外噴嘴為離心煤油噴嘴，出口直徑Df= 16.4 mm，內噴嘴為直流氧化劑噴嘴，出口直徑Do=13.4 mm?？s進室長度為10.5 mm。噴嘴安裝于模擬燃燒室中間位置，出口與模擬噴注器面板平齊。

試驗用氧化劑為400 ℃的氧氣與空氣混合形成的富氧氣，試驗用燃料為410 ℃的煤油蒸氣，采用氣氣湍流擴散燃燒可以有效模擬液氧煤油補燃循環液體火箭發動機推力室的燃燒穩定性狀態[17]。氧氣和空氣的混合物用氧化劑換熱器加熱，煤油用燃料換熱器加熱為煤油蒸氣。

1.2 試驗方法和工況

脈動壓力傳感器數目為1只，安裝在模擬燃燒室壁面距模擬噴注器面板20 mm的位置(見圖1)，型號為Kistler 6043A型，與高速數字采集系統LMS相連，采樣記錄頻率為12.8 kHz。試驗中，穩定燃料流量，逐步增加氧化劑流量至各工況點，如表1所示。混合比為氧化劑(84% 氧氣+16%空氣)與燃料(煤油)質量流量之比。

表1 試驗工況Table 1 Test operating conditions

1.3 數據處理方法

1.3.1相空間重構

相空間重構是非線性時間序列分析的基礎。對于壓力時間序列{x(i)|i=1,…,n}，由時間延遲重構技術可以得到延遲向量[20]

yi=y(i)=(x(i),x(i+τ),…,x(i+(d-1)τ)),

i=1,2,…,n-(d-1)τ

(1)

式中：d為嵌入維數，通過虛假最鄰近法(False nearest neighbor，FNN)確定；τ為延遲時間，取采樣間隔的整數倍，通過互信息法(Average mutual information，AMI)確定。

1.3.2遞歸分析

對于確定性系統，遞歸是其固有本質特性[10]。遞歸圖用來顯示和分析系統在相空間中的動力學行為。二進制矩陣R稱為遞歸矩陣，如果yi是yj的鄰近點，則矩陣元素Rij等于1，否則矩陣元素Rij等于0。數學上寫為

(2)

式中：Θ為Heaviside函數；ε為兩點間的最大距離，稱為遞歸閾值；M=n-(d-1)τ為重構相空間中向量的數量；兩條軌跡間的距離通過L2范數進行計算。遞歸圖中，白點表示0，黑點表示1。因為Rii=1，因此任何遞歸圖的主對角線上全部為黑點。

遞歸圖中線性結構的量化分析通常利用下述幾種測度。

1) 遞歸率：

(3)

遞歸率表示了黑點在遞歸圖中的密度，代表相空間中遞歸的比例分數。

2) 確定率：

(4)

確定率是遞歸圖對角結構遞歸點比例分數的測度，是系統周期行為的測度。P(l)是長度為l的對角線數量，lmin是構成對角結構的最小點數目。本文中，lmin假定為5。

3) 對角線平均長度:

(5)

平均對角線長度量化評估兩段軌跡保持在彼此鄰域的平均時間，對角線平均長度減小意味著混沌行為增加。

4) 層狀度：

(6)

層狀度表征組成豎直結構的遞歸點的比例分數。遞歸圖中，豎直結構表示層流狀態或陣發性。與確定率相似，P(v)是長度為v的豎直線數量，vmin是構成豎直結構的最小點數目。本文中，vmin假定為3。

5) 捕獲時間：

(7)

系統在某一特定狀態花費的平均時間可用捕獲時間定量描述，也是遞歸圖中豎直結構的平均長度。

上述測度在本文中用來定量表征熱聲系統的動力學特性。

2 實驗結果及分析

2.1 頻率和幅值

試驗中燃料流量不變，隨著氧化劑流量的增加(混合比增加)，系統從穩定狀態變為熱聲不穩定狀態。圖3中列出了典型工況下(混合比0.81,2.97,4.32)的壓力時間序列和FFT分析結果。從圖3可以看出，混合比0.81時系統出現了低幅值的寬頻振蕩，主要為推進劑流動激發的燃燒器低階聲腔模態，沒有發生熱聲耦合，為穩定狀態?；旌媳葹?.97和4.32時，系統為不穩定狀態，出現了多個頻率的振蕩。

圖3 不同混合比條件下的壓力時間序列和FFT分析Fig.3 Pressure time series and FFT analysis for different mixture ratios

理論計算燃燒室固有聲模態頻率時，假定各處聲速相同，縱向聲頻率按下式計算(一端開邊界、一端閉邊界)[21]

(8)

其中，c為聲速，L為燃燒室縱向長度，n為聲模態階數。切向聲頻率利用貝塞爾函數計算，其前兩階α值分別為1.8411和3.0542。計算得到模型燃燒室前五階頻率如表2所示(聲速約為1010 m/s)。試驗中出現的頻率為2736 Hz,3722 Hz,4660 Hz(混合比4.32)。由此可知，試驗中出現的不穩定頻率為燃燒室第二、第三階縱向和第一階切向聲模態頻率。根據定義[20]，它們為無公度振蕩頻率(系數為有理數時振蕩頻率線性不相關)。

表2 燃燒室前五階頻率Table 2 The first five modes of chamber

不同混合比下的最大振蕩幅值如圖4所示。在混合比為2.70時，模擬燃燒室出現了自激燃燒不穩定，混合比為2.97時，振蕩幅值達到最大，繼續增加混合比，振蕩幅值有所下降。需要注意的是，不同混合比下最大幅值對應的振蕩頻率并不相同，混合比為2.97時，第二階縱向聲模態頻率為3022 Hz。混合比變化時，火焰結構和動力學發生變化，引起局部聲阻抗和聲速的變化，因此系統振蕩頻率也會發生變化[14]。

圖4 不同混合比下最大幅值和對應的頻率Fig.4 Maximum amplitudes of FFT spectrum and corresponding frequencies at different mixture ratios

2.2 相空間重構

采用交互信息法求最優延遲時間τ時，取互信息法的第一個局部最小值為最優延遲時間[21]。求得混合比為0.81，2.97，4.32時最優延遲時間τ分別為0.11 ms,0.07 ms,0.07 ms，嵌入維數d分別為15，9，8。

利用最優延遲時間τ和最小嵌入維數d，繪制了混合比為0.81，2.97，4.32時的相空間軌跡圖(見圖5)。從圖5可以看出，混合比0.81時相軌跡圖為無序軌跡線的聚集，振蕩幅值比較低，為典型的燃燒噪聲軌跡圖?；旌媳?.97時相軌跡圖在高幅值周期振蕩(外部圓)和低幅值振蕩(內部小區域)之間連續轉變，具有陣發現象的特征[8]。混合比4.32時相軌跡圖雖然看似也為無序軌跡，但其振蕩幅值較高。由第2.1節可知，混合比4.32時其無公度振蕩頻率多于2(無公度振蕩頻率為2時相軌跡為環面[20])，因此相軌跡圖也顯示出類似無序狀態。

圖5 三種混合比下的相空間軌跡圖Fig.5 Phase space trajectory for three different mixture ratios

3 遞歸分析

3.1 遞歸圖

計算不同混合比下壓力時間序列的遞歸圖，發現混合比為0.81～2.43時的遞歸圖相近，混合比為2.7和2.97時的遞歸圖相近，混合比為3.78～4.86時的遞歸圖相近，因此下文以混合比為0.81，2.97，4.32時的遞歸圖作為典型結果進行分析。遞歸閾值ε依據吸引子直徑按比例選取，如ε=0.1，說明遞歸閾值ε取值為吸引子直徑的10%。

混合比0.81時聲壓振蕩幅值較低，為燃燒噪聲狀態，從遞歸圖6可以看出，整體呈現一種較為均勻的分布特征。進一步分析，發現存在與主對角線平行的短線，證明系統可能存在周期性過程(確定性)，而白色噪聲過程(隨機的)狀態為均勻分布的黑點[20]。

圖6 混合比k=0.81時壓力時間序列和遞歸圖Fig.6 Pressure time series and recurrence plot for mixture ratio k=0.81

混合比2.97時壓力時間曲線和對應的遞歸圖如圖7所示，壓力時間曲線上顯示高低幅值交替出現，遞歸圖上出現了密度較大和白色區域交替出現的現象。與壓力時間曲線對應可以發現，密度較大的黑色方塊對應于幅值較低的階段，白色區域對應于高幅值階段。這種白色區域和黑色方塊區域交替出現是陣發現象的典型特征。依據陣發現象的定義[23]，遞歸圖(見圖7(b))上出現的密度較大的黑色方塊對應于層流狀態(陣發之前的動力學狀態)，層流狀態被陣發狀態隔開。黑色方塊右上角上都有曲線延伸出去，形成“風箏”狀圖形，這是第II類陣發的特征[23]。出現第Ⅱ類陣發意味著系統發生了亞臨界Hopf分叉，即系統由一個穩定平衡點和一個不穩定周期解轉變為單個不穩定平衡點。發生亞臨界Hopf分叉時，系統行為與初始條件有關，存在雙穩定區或滯后區現象，在前期的試驗中得到了驗證[19]。

圖7 混合比k=2.97時壓力時間序列和遞歸圖Fig.7 Pressure time series and recurrence plot for mixture ratio k=2.97

圖8為從圖7(b)中沿對角線截取的窗口，代表著系統時間發展歷程。進入陣發之前，系統為準周期振蕩狀態b，逐漸地，系統從準周期狀態進入了周期振蕩狀態c和d。

圖8 混合比k=2.97時動力學狀態轉變過程Fig.8 Transition process of dynamic state for mixture ratio k=2.97

混合比4.32時時間曲線和遞歸圖如圖9所示。時間曲線上壓力振蕩幅值相對穩定，沒有出現混合比2.97中高低幅值振蕩交替的現象。從遞歸圖9可以看出，存在與對角線平行的分段的短線，說明系統存在周期性過程。結合第2.1節頻譜分析，此工況下系統出現了3種無公度頻率的振蕩，根據定義，系統為準周期振蕩狀態。遞歸圖符合準周期振蕩的特征：準周期振蕩會打斷較長的對角線，因此遞歸圖中會顯示出與對角線平行的分段的短線[21]。

3.2 遞歸量化分析

采用遞歸量化分析方法，計算各實驗工況下遞歸率、確定率、對角線平均長度、層狀度和捕獲時間等參數值。計算共采用時長0.78 s的數據點(10000個數據點)，分成20段(每段500個數據點)單獨進行遞歸分析，然后計算出平均值和標準差。遞歸閾值ε取為吸引子直徑的0.1[21]，計算結果如圖10所示。

遞歸率中，混合比為2.97時數值最大，對應于陣發現象，其它工況下數值相當。

圖9 混合比k=4.32時壓力時間序列和遞歸圖Fig.9 Pressure time series and recurrence plot for mixture ratio k=4.32

圖10 遞歸量化指標隨混合比的變化Fig.10 Variation of recurrence quantification analysis measures with mixture ratio

確定率中，混合比小于2.70時(燃燒噪聲狀態)確定率較低，從混合比2.70(陣發狀態)開始逐漸增加，混合比2.97時(陣發狀態)數值最大，繼續增加混合比(準周期振蕩狀態)，確定率有所降低，但仍然比燃燒噪聲階段高一個數量級。說明確定率可以將陣發狀態和準周期振蕩狀態與燃燒噪聲狀態區分出來。進一步，熱聲系統進入陣發狀態(隨后發生熱聲振蕩)后，確定率數值顯著增大，表明系統即將發生熱聲振蕩，因此，可以用確定率對熱聲振蕩進行預報。

對角線平均長度中，陣發狀態時(混合比2.70，2.97)，數值較大，說明系統的周期性較強。準周期振蕩狀態(混合比大于2.97)時，對角線平均長度與燃燒噪聲前兩個工況(混合比0.81，1.62)數值接近，原因是準周期振蕩會打斷較長的對角線[21]，減低對角線平均長度。

層狀度中，陣發狀態(混合比2.97)數值較大，驗證了層狀度對陣發現象的識別能力[21]。捕獲時間各工況下差別不大，說明各工況下局部點其動力學狀態的持續時間相當。

4 結 論

采用遞歸方法對同軸離心式噴嘴燃燒穩定性進行了分析。獲得了以下結果：

1) 試驗中燃料流量不變，增加氧化劑流量，熱聲系統依次經歷了燃燒噪聲狀態、陣發狀態和準周期振蕩狀態。

2) 燃燒噪聲狀態和準周期振蕩狀態，相軌跡圖均顯示出無序狀態；陣發狀態時，相軌跡圖在高幅值周期振蕩和低幅值振蕩之間連續轉變。

3) 燃燒噪聲狀態，遞歸圖呈現均勻分布特征。陣發狀態時，遞歸圖出現了黑色密度較大區域和白色區域交替的典型特征，進一步判斷為第II類陣發現象，意味著系統發生了亞臨界Hopf分叉，與前期實驗結果一致。準周期振蕩狀態，遞歸圖顯示出與對角線平行的分段的短線。

4) 遞歸量化分析表明，確定率可以將陣發狀態和準周期振蕩狀態與燃燒噪聲狀態區分出來，并可以對熱聲振蕩進行預報。