基于變分模態(tài)分解的燃燒室熱聲振蕩衰減系數(shù)辨識(shí)

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

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

0 引言

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

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

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

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

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

1.1 燃燒室熱聲模型

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

(1)

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

(2)

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

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

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

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

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

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

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

(3)

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

(4)

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

(5)

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

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

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

(6)

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

(7)

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

(8)

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

(9)

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

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

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

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

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

(10)

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

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

(11)

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

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

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

(12)

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

(13)

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

(14)

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

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 燃燒室動(dòng)態(tài)壓力監(jiān)測(cè)圖

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

圖6 燃燒室壓力信號(hào)功率譜密度

圖7為VMD算法對(duì)脈動(dòng)壓力數(shù)據(jù)的處理結(jié)果,經(jīng)過計(jì)算:分離模態(tài)數(shù)K取值為9,噪聲容限τ為0.92,共分離出9個(gè)本征模態(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)相近,且與壓力信號(hào)功率譜密度圖中2個(gè)寬頻凸峰相近(2.534 kHz、4.814 kHz),可以看出VMD算法可準(zhǔn)確地分離得到諧振模態(tài)中心頻率。

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

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

表3 VMD辨識(shí)結(jié)果

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

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

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

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

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

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

3 結(jié)論

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

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

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