楊尚榮, 吳林龍, 于 涵, 楊寶娥
(西安航天動力研究所 液體火箭發動機技術重點實驗室, 西安 710100)
燃燒不穩定是航空航天發動機研制中需要重點關注的問題,通常表現為燃燒室聲模態頻率附近的高幅值壓力振蕩。國內外在其產生機理、評估和預測方法、以及控制手段等方面開展了大量的研究[1-2]。目前預測方法在工程上的應用還未成熟,尚不能在發動機試驗前對其穩定性進行可靠的預測。因此一旦出現燃燒不穩定,通常從兩方面出發去解決:一方面減少激勵能量;另一方面增加燃燒室聲學阻尼。前者需要考慮燃燒和聲場間的耦合作用過程,由于燃燒現象本身的復雜性,實施改進方案難度較大。相比較而言,被動阻尼裝置的設計要容易一些。
從線性角度,增長率表征小擾動下燃燒室振蕩能量的變化率。如果其值大于零,說明系統不穩定,反之,說明系統穩定,此時增長率也稱為衰減率。對于出現燃燒不穩定的燃燒室,新設計阻尼裝置的衰減率應大于燃燒室的增長率并保留一定的裕度,即需要燃燒室不穩定聲模態的增長率作為設計輸入參數。
學者們提出了多種方法來識別增長率。主要分為外部激勵方法(輸入-輸出識別)和無外部激勵方法(僅輸出識別)。外部激勵方法中,對于Hopf類型的分叉,Lee等[3]利用分叉前的噪聲相干共振特性[4],識別系統的非線性燃燒響應,外推出系統分叉點的位置和類型、增長率、極限環振蕩幅值等,進一步可以對系統穩定性邊界進行預測[5]。另一類方法采用主動控制[6],在不穩定工況點施加和關閉控制,獲得振蕩幅值從小擾動增加到極限環的過程。利用振蕩幅值初始增長階段(滿足線性假設)數據,擬合指數函數或采用其他數據處理方法(如DMD(dynamic mode decomposition)[7])識別增長率。實際發動機試驗中加入特定形式的外部激勵或主動控制比較困難,因此外部激勵方法在工程應用中不方便。
無外部激勵的方法,第一類與主動控制方法中關閉控制后的處理方法相同,區別是利用自發燃燒不穩定振蕩幅值初始增長階段數據[8]。另一類方法基于噪聲激勵作用下燃燒室聲模態的隨機動力學理論[9],從極限環振蕩數據中提取線性增長率。根據噪聲強度等級,Noiray等[10]提出對應的3種識別方法,并利用燃氣輪機中壓力振蕩數據進行了比較。進一步,利用數值模擬[11]和燃燒器試驗[12]對第三種方法進行了驗證。Bonciolini等[13]比較了彩色噪聲和白噪聲的激勵效果,發現在振蕩頻率附近采用合適的帶通濾波,兩者的區別可以忽略。若帶通濾波過窄,則兩者的識別結果差別較大。Boujo等[14]通過引入伴隨優化方法解決了上述問題。即使考慮燃燒響應相對于壓力振蕩的時間延遲,該方法依然可以較準確的識別線性增長率[15]。Hummel等[16]將方法推廣到周向不穩定對聲模態增長率的識別。
本文開展同軸離心式噴嘴穩定性試驗研究,采用無外部激勵方法識別不同工況下的線性增長率,并對其參數敏感性開展分析,為液體火箭發動機燃燒穩定性分析提供參考。
燃燒試驗系統的詳細介紹參見文獻[17]。模擬燃燒室為敞口圓筒形結構。試驗用噴嘴(如圖1所示)為帶縮進室的同軸離心噴嘴,外噴嘴為煤油離心噴嘴,內噴嘴為直流氧化劑噴嘴。

1. 直流噴嘴; 2. 離心噴嘴; 3. 節流嘴; 4. 富氧氣通道; 5. 縮進室。圖1 帶縮進室的同軸離心噴嘴Fig.1 Recessed swirl coaxial injector

表1 試驗工況Tab.1 Test operating conditions
當燃燒室只有單階聲模態振蕩時,可以不考慮模態間的耦合作用。壓力振蕩可以用非線性振動方程描述
式中:η為脈動壓力;f為由氣動力和燃燒過程產生的激勵源項;α為聲模態衰減率;ω0為聲模態角頻率;ξ為燃燒噪聲,假定為高斯白噪聲,有〈ξξτ〉=Γδ(τ),Γ為噪聲強度。

(2)
假定壓力振蕩接近諧振,可以將脈動壓力寫成如下幅頻形式
η(t)=A(t)cos[ωt+σ(t)]
(3)
式中,A(t)和σ(t)分別為壓力振蕩的幅值和相位。代入式(2),利用隨機平均方法,式(2)可以寫成幅值A和相位σ的一階隨機微分方程。由于A和σ可以解耦,僅列出幅值A滿足的方程
(4)

(5)
駐定情況下,FP方程存在解析解
(6)



(7)
進一步通過非線性擬合識別出線性增長率ν、飽和系數κ和噪聲強度Γ。
不同混合比條件下的壓力時間序列及其對應的相空間軌跡如圖2所示。當混合比為15.0時,壓力振蕩幅值較小,認為是燃燒噪聲階段,相軌跡圖為無序軌跡線的聚集。當混合比為17.5時,振蕩幅值增加,振幅波動也較大,軌跡在不同幅值周期振蕩之間轉變,中心為空心。當混合比為20.0和22.5時,達到極限環振蕩,軌跡圖為圓環面,環的寬帶代表噪聲對幅值的影響。時頻分析顯示振蕩燃燒頻率在(2 340±100)Hz內。

圖2 壓力時間序列和相空間軌跡圖Fig.2 Pressure time series and phase portrait
試驗中壓力數據有單位,如Pa,kPa等。本節分析壓力數據取不同單位時對識別結果的影響。令幅值A(t)=kB(t),k為單位間的換算關系,B(t)為換算后的幅值時間序列。代入式(5),化簡并結合可得
(8)
比較式(5)和式(8),可知單位變換后,增長率不變,飽和系數為原來的k2倍,噪聲強度為原來的1/k2倍。

表2 不同幅值下參數識別結果
本文壓力單位選用kPa進行數據分析,其他單位下的飽和系數和噪聲強度可通過2.2節換算關系得到。4種工況下壓力時間序列識別結果如圖3所示。隨著混合比增加(總流量也增加),壓力振蕩均方根(root mean square,RMS)幅值增加,增長率ν由負值轉變為正值,飽和系數κ下降,噪聲強度Γ增加。具體地,當混合比為15.0時,其增長率為負值,系統是穩定的。當混合比為17.5時,增長率也為負值,從線性角度分析系統應該是穩定的。但從圖2壓力曲線和相軌跡可知,其周期性振蕩幅值也較大。說明存在噪聲條件下(實際發動機中燃燒噪聲不可避免),不僅要保證系統增長率為負,還需要留有一定裕度,才能保證系統不發生大幅值周期振蕩。當混合比為20.0和22.5時,增長率均為正值,且幅值越高,增長率越大。該結果不能作為通用結論來推廣,理論上,以式(2)作為控制方程的系統,壓力振蕩幅值由增長率ν和飽和系數κ共同決定。不同工況下,飽和系數κ非定值。此外試驗中也發現較高幅值的壓力振蕩并非一定對應較大的增長率[20]。

圖3 試驗壓力時間序列RMS值和識別結果Fig.3 RMS values of test pressure time series and identification results
將識別結果代入FP方程解析解式(6),與試驗壓力振蕩幅值概率密度分布進行比較,如圖4所示。當混合比為20時,試驗結果和識別結果差別較大,其他工況結果相對較好,一定程度上驗證了燃燒響應模型和識別方法的合理性。Oriray研究中的數值試驗發現增長率在-5~5 rad/s內時,準確識別需要約30 s時長的數據。本文試驗采樣時間T=5 s,相比準確識別所需的數據長度而言較短,因此需要對識別結果的不確定性進行分析。

圖4 振蕩幅值概率密度函數Fig.4 Probability density function of oscillation amplitude
本節分析識別結果的不確定性以及數據采樣時間和采樣頻率對識別結果的影響。每個工況下試驗數據只有一次,且系統真實的增長率未知。因此,為了量化識別方法的不確定性,將每個工況下的識別值作為設定值,數值求解隨機微分式(2)生成100組壓力時間序列,識別其增長率并取平均值和標準差,來量化該工況下識別結果的不確定性。
4個工況下的計算結果如圖5所示。采樣時間T和采樣率sf與試驗中取值相同,即T=5 s,sf=25.6 kHz。結果表明,當混合比為20.0和22.5時,增長率均值與設定值偏差分別為11%和15%,標準差分別為17%和20%。當混合比為15.0和17.5時,增長率均值與設定值偏差分別為76%和50%,標準差分別為62%和56%。工程應用時可以確定一個閾值,設定值(試驗數據識別值)與識別值偏差小于該閥值時,則可以將識別值的標準差作為設定值的標準差,從而確定識別結果的不確定性。

圖5 采樣率對增長率識別的影響,采樣時間T=5 sFig.5 Effect of sampling rate on growth rate identification, sampling time T=5 s
采樣頻率對結果的影響見圖5,采樣時間均為5 s,采樣率分別為25.6 kHz,50 kHz,1 000 kHz。均值和標準差均由式(2)生成的100組壓力時間序列計算得到。對識別值的橫坐標作了平移,以便圖形能清晰表達。增加采樣率,當混合比為15.0和17.5的偏差(設定值與識別均值之差)減小,混合比20.0和22.5的偏差在采樣率50 kHz時最小,4個混合比下其識別值的標準差均增加。
采樣時間對結果的影響如圖6所示,采樣頻率均為1 000 kHz,采樣時間分別為1 s,3 s,5 s。每個均值和標準差均由式(2)生成的100組壓力時間序列計算得到。同樣對識別值的橫坐標作了平移。從圖6可知,增加采樣時間,4個混合比下識別值的偏差和標準差均減小。由此可知,為了增加識別的精度,應該增加采樣時間,同時提高采樣頻率,以便在數據分析時可選出對應工況下較優的采樣頻率。

圖6 采樣時間對增長率識別的影響,采樣率sf=1 000 kHzFig.6 Effect of sampling time on growth rate identification, sampling rate sf=1 000 kHz
增長率也可以從壓力振蕩幅值線性增長階段識別得到。由于壓力振蕩幅值在線性增長階段以指數形式增長[21],因此可對壓力振蕩時間序列取包絡,利用指數函數去擬合,或將包絡取對數后利用線性函數去擬合。本章分析噪聲對增長率識別結果的影響。


圖7 無噪聲壓力時間序列和增長率識別Fig.7 Noise free pressure time series and growth rate identification

圖8 含噪聲壓力時間序列和增長率識別Fig.8 Pressure time series with noise and growth rate identification

圖9 噪聲對增長率識別的影響Fig.9 Effect of noise on growth rate identification
(1)壓力振蕩幅值的單位不改變增長率的值,但會影響飽和系數和噪聲強度,理論上得到了相似關系并獲得了數值試驗的驗證。
(2)即使增長率小于零,如果噪聲強度較大,系統也有可能發生較大幅值的壓力振蕩,因此系統增長率應該留有一定裕度。
(3)提高采樣時間可以提高識別精度,不同工況下存在各自較優的采樣率。
(4)當噪聲強度較大時,從振幅線性增長階段識別增長率存在較大誤差,需要考慮噪聲導致的不確定性。