樊志偉,鐘 戰(zhàn),李 庚,聶萬勝,何 博

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

0 引言

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

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

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

本文基于隨機白噪聲驅(qū)動的二階非線性振子模型,對理論模型仿真數(shù)據(jù)和針?biāo)ㄊ侥Ｐ桶l(fā)動機熱試車數(shù)據(jù)進行分析,重點比較了現(xiàn)有帶通濾波方法與變分模態(tài)分解算法的優(yōu)劣,驗證了變分模態(tài)分解算法的適用性與準(zhǔn)確性,并將其應(yīng)用于模型發(fā)動機試驗結(jié)果處理,為液體火箭發(fā)動機燃燒穩(wěn)定性裕度評估提供參考。

1 系統(tǒng)模型與辨識方法

1.1 燃燒室熱聲模型

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

(1)

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

(2)

式(2)中,當(dāng)v<0時系統(tǒng)處于線性穩(wěn)定狀態(tài)。由此可知,根據(jù)辨識得出的衰減系數(shù)v與0之間的距離即可評定燃燒裝置的燃燒穩(wěn)定性裕度。

1.2 模態(tài)提取方法

衰減系數(shù)辨識方法中最重要的就是特征模態(tài)的獲取,而前期研究表明:帶通濾波參數(shù)選取對辨識精度有較大影響,特別是當(dāng)處于穩(wěn)定燃燒階段,無法精準(zhǔn)獲取燃燒室內(nèi)熱聲耦合頻率,進而影響濾波寬度的選擇[13]。現(xiàn)有辨識方法采用的是將燃燒室固有聲學(xué)頻率作為特征模態(tài)中心頻率,進而測算合適的濾波寬度,以此來獲取特征模態(tài)信息。

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

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

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

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

(3)

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

(4)

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

(5)

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

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

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

(6)

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

(7)

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

(8)

5)重復(fù)循環(huán)步驟2)～4),直到滿足如下約束條件

(9)

式中:τ為噪聲容限;K為預(yù)設(shè)分解層數(shù)。

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

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

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

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

(10)

1.3 時/頻域衰減系數(shù)辨識方法

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

(11)

式中:ωi為特征模態(tài)的角頻率。

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

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

(12)

則可得到幅值方程的自相關(guān)函數(shù)為

(13)

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

(14)

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

2 結(jié)果與討論

2.1 模態(tài)提取與算法驗證

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

表1 數(shù)值模擬信號系統(tǒng)參數(shù)

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

圖1 數(shù)值模擬信號及其頻譜圖

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

圖2 中心頻率對系統(tǒng)辨識結(jié)果影響

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

圖3 濾波寬度對系統(tǒng)辨識結(jié)果影響

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

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

圖4 VMD分離模態(tài)結(jié)果

表2 VMD與帶通濾波方法辨識結(jié)果

2.2 試車數(shù)據(jù)與驗證

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

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

圖5 燃燒室動態(tài)壓力監(jiān)測圖

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

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

圖7為VMD算法對脈動壓力數(shù)據(jù)的處理結(jié)果,經(jīng)過計算:分離模態(tài)數(shù)K取值為9,噪聲容限τ為0.92,共分離出9個本征模態(tài)函數(shù),中心頻率分別為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,其中模態(tài)2、模態(tài)4的中心頻率(2.510 kHz、4.882 kHz)與燃燒室固有2階縱向聲頻率(2.528 kHz)和4階縱向聲頻率(5.056 kHz)相近,且與壓力信號功率譜密度圖中2個寬頻凸峰相近(2.534 kHz、4.814 kHz),可以看出VMD算法可準(zhǔn)確地分離得到諧振模態(tài)中心頻率。

圖7 VMD分離模態(tài)結(jié)果

進而對IMF2和IMF4進行衰減系數(shù)辨識,辨識結(jié)果如表3所示。

表3 VMD辨識結(jié)果

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

圖8 濾波寬度對系統(tǒng)辨識結(jié)果影響

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

表4 VMD與帶通濾波方法辨識結(jié)果

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

圖9 衰減系數(shù)隨時間變化

3 結(jié)論

本文開展了基于變分模態(tài)分解的燃燒室熱聲振蕩衰減系數(shù)辨識研究,基于理論模型仿真數(shù)據(jù)進行驗證,并應(yīng)用在模型發(fā)動機熱試車數(shù)據(jù)處理中,得出以下結(jié)論。

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

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