基于變分模態分解的燃燒室熱聲振蕩衰減系數辨識

樊志偉,鐘 戰,李 庚,聶萬勝,何 博

(航天工程大學 宇航科學與技術系,北京 101400)

0 引言

高頻燃燒不穩定源自燃燒室內燃燒過程與聲學振蕩的耦合[1-2],一旦發生將在幾秒內破壞整個發動機。為盡量避免液體火箭發動機使用過程中出現燃燒不穩定,提高工作可靠性,通常需要在其研制階段開展大量燃燒穩定性裕度試驗評估工作。目前,評估燃燒穩定性裕度主要基于如下兩種技術指標參數[3-4]:①燃燒室工作過程對人為擾動響應的馳豫時間[5];②燃燒室容腔諧振頻率的振蕩衰減率。其中,前者是國內外早期試驗評估通常使用的方法,但需要外部激勵裝置誘導激勵燃燒室出現高幅值壓力振蕩,會對發動機結構產生不利影響;后者則僅需要發動機穩定燃燒階段的燃燒室壓力監測數據,是當前極具發展前景的一種燃燒穩定性裕度評估方法。

對于燃燒室熱聲振蕩衰減率的辨識,其基本原理是將其穩定燃燒階段的系統動力學建模為不同復雜程度振子模型方程[6],進而通過系統辨識方法從燃燒室脈動壓力數據中辨識出模型方程衰減系數,從而判定發動機燃燒穩定性裕度。Lieuwen利用壓力數據的自相關函數計算出燃燒系統衰減系數,從時域對非線性耦合振子模型的衰減系數進行了辨識[7]。Stadlmair等采用貝葉斯統計方法,從單特征模態拓展到了多模態的衰減系數辨識[8]。Yi等利用維納辛欽定理將時域方法變換到了頻域,擬合脈動壓力功率譜得到衰減系數[9]。Noiray等基于噪聲相干共振現象,將系統模型建模為隨機噪聲激勵的范德波爾方程,利用壓力信號概率密度函數建立了系統模型的?？?普朗克方程,最終通過求解方程的Kramers-Moyal系數擬合出系統衰減系數[10-11]。Boujo等通過引入伴隨優化方法修正了Noiray方法中因有限采樣時間效應帶來的辨識誤差,進一步提高了辨識精度[12]。

然而,衰減系數辨識方法中誤差影響最大的環節是獲取燃燒室熱聲諧振模態。目前對于燃燒室內諧振模態信息的提取主要采用帶通濾波方法,這種方法需要預先確定帶通濾波中心頻率與濾波寬度。但是,燃燒室穩定噪聲燃燒階段的中心頻率通常難以通過理論模型獲取,因此使用理論聲學諧振頻率近似替代中心頻率會存在較大誤差,從而影響辨識結果精度。此外,帶通濾波寬度的選取也會對衰減系數辨識結果產生重要影響[13]。楊尚榮等提出了一種基于非線性最小二乘法的濾波參數計算方法,然而該方法僅適用于單個諧振模態的分離辨識,對兩鄰近模態混疊情況下分離辨識仍存在不足[14]。

本文基于隨機白噪聲驅動的二階非線性振子模型,對理論模型仿真數據和針栓式模型發動機熱試車數據進行分析,重點比較了現有帶通濾波方法與變分模態分解算法的優劣,驗證了變分模態分解算法的適用性與準確性,并將其應用于模型發動機試驗結果處理,為液體火箭發動機燃燒穩定性裕度評估提供參考。

1 系統模型與辨識方法

1.1 燃燒室熱聲模型

目前對燃燒穩定性研究表明:燃燒不穩定是火焰、流場和聲場間非線性耦合的結果,而液體火箭發動機的高頻燃燒不穩定大多是火焰與聲場之間的熱聲耦合導致的[1]。所以,燃燒不穩定低階動力學模型通常主要考慮聲場與火焰熱釋放之間的耦合,而將湍流等流場影響簡化為隨機白噪聲項,則熱聲耦合振蕩模型可用如下隨機振子方程描述。

(1)

式中:p′(x,t)為時間t位置x處的脈動壓力;ηi和ψi分別為第i階聲模態的幅值和振型;Dε(t)為燃燒室壓力振蕩中湍流引起的白噪聲擾動,D為噪聲強度,ε(t)為單位高斯白噪聲;f為火焰熱釋放對壓力振蕩幅值的反饋函數;αi為聲模態衰減系數;ω為角頻率。

(2)

式(2)中,當v<0時系統處于線性穩定狀態。由此可知,根據辨識得出的衰減系數v與0之間的距離即可評定燃燒裝置的燃燒穩定性裕度。

1.2 模態提取方法

衰減系數辨識方法中最重要的就是特征模態的獲取,而前期研究表明:帶通濾波參數選取對辨識精度有較大影響,特別是當處于穩定燃燒階段,無法精準獲取燃燒室內熱聲耦合頻率,進而影響濾波寬度的選擇[13]?，F有辨識方法采用的是將燃燒室固有聲學頻率作為特征模態中心頻率,進而測算合適的濾波寬度,以此來獲取特征模態信息。

然而燃燒噪聲階段,燃燒室壓力監測信號為寬頻振蕩,較難選擇中心頻率。另外帶通濾波寬度也對衰減系數辨識有嚴重影響。當帶通濾波寬度較小時,帶通濾波提取的模態信息會丟失部分熱聲耦合模態從而導致辨識精度下降;而帶通濾波寬度過大時,帶通濾波提取的模態信息會包含多個特征模態信息從而導致辨識精度下降。所以得到準確可靠的熱聲耦合模態信息不僅要求包含確定的中心頻率,還要求提取模態信息在頻譜上存在一定寬度,以包含所有與該模態耦合的湍流噪聲信息。

隨著信號處理方法的發展,經驗模態分解[15](empirical mode decomposition,EMD)、變分模態分解[16](variational mode decomposition,VMD)、局部均值分解[17](local mean decomposition,LMD)等方法逐漸應用于時間序列分解中,然而EMD、LMD方法都存在模態混疊、邊界效應和噪聲魯棒性弱等不足,而VMD方法雖然具有良好的噪聲魯棒性,但由于其需要預先設定好分解層數,致使分解結果可能存在過分解或欠分解情況。本文采用文獻[18]提出的基于峭度的VMD分解方法計算分解模態中的最大相關分量的峭度值,當峭度值達到最大時停止分解,最終得到最優模態分解結果。

變分模態分解方法將原始信號分解成n個限制帶寬的本征模態函數(intrinsic mode function,IMF),表示為

un(t)=An(t)cos[ωn(t)]

(3)

式中:un(t)為分解得到的n個IMF;An(t)為un(t)的瞬時幅值;ωn(t)為un(t)的瞬時頻率。每個分量在中心頻率處集中,可用高斯平滑調制信號來估算帶寬,由于VMD分解的稀疏性,分解時可將其轉換為

(4)

為解決上述的約束最優化問題,將約束變分問題轉換為非約束變分問題,添加二次懲罰項和拉格朗日因子,得到如下增廣拉格朗日函數以求上述問題的最優解。

(5)

式中:α為乘法因子;λ為拉格朗日因子。

利用變分模態分解方法將原始信號分解為K個IMF,過程如下。

1)初始化{uk}、{ωk}、λ1為0。

(6)

3)對所有ω≥0,更新泛函ωk,即

(7)

4)對所有ω≥0,更新拉格朗日因子λ,即

(8)

5)重復循環步驟2)～4),直到滿足如下約束條件

(9)

式中:τ為噪聲容限;K為預設分解層數。

從分解過程可以看出分解層數與噪聲容限對分解結果影響最大?；谇投鹊姆纸鈱訑礙確定方法是以最大相關IMF的峭度作為對照依據,以峭度最大作為循環截止條件來確定K的最優解。具體步驟如下。

2)計算ui(t)與原始信號x(t)的相關系數,得到最大相關分量uj(t)。

4)重復步驟1)～3),直到找到峭度最大時的分解層數K。

在確定分解層數后,利用殘差指數REI確定噪聲容限,在τ∈[0,1]范圍下選擇殘差最小的噪聲容限τ。殘差指數計算式為

(10)

1.3 時/頻域衰減系數辨識方法

對式(2)進行傅里葉變換并取模的平方后,可得到幅值方程的功率譜函數為

(11)

式中:ωi為特征模態的角頻率。

系統參數D、v可由壓力時間序列辨識得到,即頻域衰減系數辨識方法。

根據維納—辛欽定理可得時域計算的自相關函數與頻域計算的功率譜Sη(ω)的關系為

(12)

則可得到幅值方程的自相關函數為

(13)

由于線性衰減系數v遠小于角頻率ωi,則自相關函數第二項可忽略,由此可得

(14)

由式(14)可知自相關函數幅值在頻率ωi附近振蕩,且幅值成指數衰減,衰減率為線性衰減系數v,由此可得線性衰減系數,即時域的衰減系數辨識方法。

2 結果與討論

2.1 模態提取與算法驗證

理論驗證信號采用非線性耦合振子模型,模型方程為式(2),設定模型參數如表1所示,需辨識的信號由3種不同衰減率模擬信號疊加而成,數值模擬時間為10 s,采樣頻率為10 kHz。

表1 數值模擬信號系統參數

圖1為3種指定參數模擬信號疊加后的幅值與頻譜圖。

圖1 數值模擬信號及其頻譜圖

圖2為帶通濾波中心頻率選擇影響,當帶通濾波寬度選取為中心頻率的20%時,改變帶通濾波中心頻率,可見隨著中心頻率選取從1.0～1.2 kHz,頻域方法辨識結果越來越向模態3靠近;時域方法辨識結果先減小,然后趨于模態3的預設衰減系數值(-0.05),誤差越來越大,逐漸偏離相對應的模擬信號預設衰減系數。其中,時域方法先降低是因為2個相近模態(0.9 kHz、1.0 kHz)之間的混疊效應。

圖2 中心頻率對系統辨識結果影響

圖3為帶通濾波寬度對模擬信號的辨識結果影響,可見時域方法辨識結果先增大后減小,最后仍存在一定誤差,頻域方法辨識結果一直比較穩定。這是因為時域方法受濾波寬度影響較大,起初濾波寬度較小,模態2(1.0 kHz)只有部分信息進行辨識;進而因模態1(0.9 kHz)與模態2(1.0 kHz)兩者中心頻率相近,濾波寬度增大會將模態1的影響引入辨識結果中;最后隨著濾波寬度繼續加大,模態3(1.2 kHz)也影響著辨識結果。頻域方法辨識結果較好是因為中心頻率不變,濾波寬度影響程度較小。

圖3 濾波寬度對系統辨識結果影響

在傳統模態提取方法中,頻域方法對中心頻率的選取比較敏感,時域方法則需要合適的中心頻率與濾波寬度。然而在燃燒噪聲階段,燃燒室壓力表現為以固有頻率為中心的寬頻振蕩,燃燒噪聲自身的隨機性會對辨識結果有較大影響。因此,準確可靠的模態提取方法不僅需要準確的中心頻率,也需要一定的頻帶寬度以包含所有模態信息。

圖4為本文VMD算法處理的結果。處理VMD分解得到的結果如表2所示,可見VMD算法可以高效分離得到中心頻率與模態信號,所得的模態經過時/頻域方法辨識后與帶通濾波結果誤差較小,可以使用基于峭度最大的變分模態分解算法處理壓力信號,與帶通濾波方法相比,變分模態分解可以自動計算得到中心頻率及信號,無需先驗的熱聲諧振頻率。

圖4 VMD分離模態結果

表2 VMD與帶通濾波方法辨識結果

2.2 試車數據與驗證

試驗發動機采用氣氧/乙醇作為推進劑,點火過程利用火花塞來啟動發動機,噴嘴類型為針栓式,試驗工況設置為:氧化劑(O2,300 K)為112 g/s,燃料(質量分數為50%的乙醇水溶液,300 K)為218 g/s,氧燃比為0.51;該工況下燃燒室內總溫1 506 K,聲速為910 m/s,模型燃燒室為圓筒形結構,燃燒室直徑為70 mm,長度為350 mm,收縮段長度為20 mm,通過計算可得1階縱向聲學固有頻率為1.264 kHz,2階 為2.528 kHz,4階為5.056 kHz。壓力傳感器選用一個采樣頻率為125 kHz的Kistler 6043A動壓傳感器,脈動壓力測點距噴注面板平面15 mm。

圖5為燃燒室脈動壓力監測信號圖,從6.1 s開始進行點火,在6.2～7.2 s之間為燃燒階段,7.2 s時刻停止供應氧氣,燃燒持續時間約1.0 s,期間燃燒室內壓力約為1.0 MPa,由放大圖可知,發動機處于穩定燃燒狀態。為保證實驗安全,發動機點火在富氧條件下進行,熄滅在富燃條件進行。

圖5 燃燒室動態壓力監測圖

圖6為燃燒室壓力功率譜密度,可見在穩定燃燒階段(6.2～7.2 s)存在2個凸頻,在頻率2.534 kHz、4.814 kHz附近顯示寬頻燃燒振蕩,分別對應理論計算出的燃燒室2階縱向、4階縱向聲模態。

圖6 燃燒室壓力信號功率譜密度

圖7為VMD算法對脈動壓力數據的處理結果,經過計算:分離模態數K取值為9,噪聲容限τ為0.92,共分離出9個本征模態函數,中心頻率分別為0.195 kHz、2.510 kHz、3.710 kHz、4.882 kHz、7.421 kHz、10.156 kHz、20.210 kHz、22.613 kHz、32.738 kHz,其中模態2、模態4的中心頻率(2.510 kHz、4.882 kHz)與燃燒室固有2階縱向聲頻率(2.528 kHz)和4階縱向聲頻率(5.056 kHz)相近,且與壓力信號功率譜密度圖中2個寬頻凸峰相近(2.534 kHz、4.814 kHz),可以看出VMD算法可準確地分離得到諧振模態中心頻率。

圖7 VMD分離模態結果

進而對IMF2和IMF4進行衰減系數辨識,辨識結果如表3所示。

表3 VMD辨識結果

圖8為濾波寬度對線性衰減系數影響,隨著濾波寬度增加,衰減系數隨之減小。所以在帶通濾波寬度為中心頻率的30%時,衰減系數變化基本穩定,因此時域、頻域方法中帶通濾波寬度均選擇為中心頻率的30%。

圖8 濾波寬度對系統辨識結果影響

進而對濾波結果進行衰減系數辨識,所得結果與VMD方法所得結果對比如表4所示,對比結果可得無論2階還是4階模態的衰減系數,時域方法所得結果相差較小,頻域方法所得結果差異相對較大,具體表現在基于VMD方法濾波后的頻域結果更高,主要是因為VMD濾波所得的頻帶寬度更大。

表4 VMD與帶通濾波方法辨識結果

圖9為衰減系數隨時間變化圖,對發動機穩定燃燒階段數據進行滑窗切片處理(滑窗寬度為0.20 s,前進步長為0.01 s),對切片數據進行衰減系數辨識,由此可見衰減系數在穩定燃燒階段保持穩定,2階縱向模態在開機表現出較強的不穩定趨勢,4階縱向模態在啟停期間并未有明顯不穩定趨勢。出現這種現象的原因是在啟動點火階段,燃燒室燃料分布不均勻,還未形成穩定燃燒,所以未能形成熱聲耦合共振場[19-20],得到的衰減系數本質上反映了燃燒室結構的頻率選擇性。

圖9 衰減系數隨時間變化

3 結論

本文開展了基于變分模態分解的燃燒室熱聲振蕩衰減系數辨識研究,基于理論模型仿真數據進行驗證,并應用在模型發動機熱試車數據處理中,得出以下結論。

1)基于最大峭度準則的變分模態分解可以減小現有辨識方法中中心頻率選擇的任意性,有效得出燃燒室熱聲諧振頻率,提高計算精度。

2)模型發動機試車工況為穩定燃燒工況,2階縱向模態衰減系數為-0.35,4階縱向模態衰減系數為-0.13,同時變分模態分解算法與帶通濾波方法所得結果差異體現在頻域方法結果,是由變分模態分解算法濾波后頻帶寬度較大所致。