趙 瑜，李莎莎，劉 喆，薛牧遙

(上海航天動力技術研究所，上海 201109)

0 引言

在發(fā)動機工作過程中，燃燒室受到擾動時就會產(chǎn)生一定振型和頻率的聲波。當擾動產(chǎn)生的壓力振蕩頻率與燃燒室聲場的固有頻率一致時，將產(chǎn)生共振現(xiàn)象，即聲不穩(wěn)定燃燒。此時，燃燒室內(nèi)壓強的振蕩幅值將明顯放大，從而導致發(fā)動機產(chǎn)生劇烈的機械振動、旋轉，并破壞燃燒室殼體，甚至可能引起發(fā)動機爆炸并形成災難性后果。

目前，固體火箭發(fā)動機裝藥形式越來越復雜，很多發(fā)動機采用分段式裝藥，段與段之間存在絕熱環(huán)等障礙物，使燃氣在流動過程中發(fā)生分離并出現(xiàn)漩渦，而固體發(fā)動機燃燒室內(nèi)規(guī)律性的漩渦脫落被認為是引起聲渦耦合的重要原因。聲渦耦合理論指出：當流場中漩渦脫落頻率接近燃燒室聲腔的某一階聲模態(tài)的固有頻率時，就有可能激發(fā)該模態(tài)并引起聲場共振，導致燃燒室內(nèi)壓強振蕩，出現(xiàn)不穩(wěn)定燃燒現(xiàn)象[1]。國內(nèi)外眾多學者對此現(xiàn)象展開過研究。如：JACOB等[2]首次提出了聲渦耦合是引起固體火箭發(fā)動機不穩(wěn)定燃燒的潛在誘因；DOSTON等[3]研究了Tiatan IV RSRM的推力振蕩特性，建立漩渦形成、脫落、碰撞、聲反饋的4步模型，為固體火箭發(fā)動機渦聲耦合研究奠定了基礎；CULICK等[4]通過試驗證實了旋渦脫落頻率與發(fā)動機聲腔內(nèi)的聲頻率可相互耦合，當兩者頻率接近時會出現(xiàn)壓力振蕩。文獻[5-7]用P230縮比冷流試驗發(fā)動機對障礙物漩渦脫落進行了大量試驗研究，并開展了相應的數(shù)值模擬工作，對渦聲耦合引起的壓力振蕩現(xiàn)象進行了深入分析。國內(nèi)以西北工業(yè)大學和北京理工大學為主，在不穩(wěn)定燃燒的理論研究、試驗測試、數(shù)值模擬等方面開展了一系列研究。文獻[8-11]相繼對固體火箭發(fā)動機內(nèi)的聲渦耦合現(xiàn)象進行了數(shù)值模擬，研究表明：大渦模擬技術能較好地捕捉燃燒室內(nèi)漩渦的運動規(guī)律及壓力振蕩特性，是不穩(wěn)定燃燒研究的有力工具。文獻[12-14]詳細分析了聲渦耦合壓力振蕩現(xiàn)象，研究了固體發(fā)動機工作末期不穩(wěn)定燃燒的機理，并提出了頭部空腔抑制壓力振蕩的方法。蘇萬興等[15]以VKI發(fā)動機為模型，采用大渦模擬開展了障礙物旋渦脫落引起的壓力振蕩現(xiàn)象研究，通過數(shù)值模擬、試驗研究與理論分析相結合，系統(tǒng)研究了大長徑比固體發(fā)動機工作穩(wěn)定性的關鍵增益機理與阻尼特性。

本文以某雙脈沖發(fā)動機為研究對象，探尋其二脈沖工作時壓力振蕩大幅上升的原因，以及擾流環(huán)抑振效果的內(nèi)在機理。經(jīng)初步分析，發(fā)動機燃燒室聲腔的一階軸向固有頻率(由經(jīng)典圓柱形聲腔固有頻率預估公式得到)與其內(nèi)部的壓力振蕩頻率(由傅里葉變化得到)十分接近，初步判斷可能是聲渦耦合誘發(fā)的不穩(wěn)定燃燒現(xiàn)象。為驗證該推測，采用有限元和大渦模擬算法分別對燃燒室聲腔固有頻率和漩渦脫落頻率進行了深入分析，并在此基礎上研究了使用擾流環(huán)后帶來的影響。

1 計算模型與數(shù)值算法

1.1 有限元聲振型計算方法

固體火箭發(fā)動機燃燒室空腔內(nèi)充滿燃氣，受到初始擾動時會出現(xiàn)各種振型的自由聲振蕩，在聲波的傳播過程中，介質(zhì)質(zhì)點的振動會導致局部壓力和密度在穩(wěn)定值附近振蕩，可表示為

(1)

(2)

式中：a為聲速。

有限元模型通用的矩陣方程為

(3)

式中：M為結構質(zhì)量矩陣；C為結構阻尼矩陣；K為結構剛度矩陣。為獲得燃燒室空腔的振型及相應的頻率，需要根據(jù)聲壓與聲速關系，將式(2)轉換為如下單元矩陣形式

(Kf-ω2Mf)p=0

(4)

式中：Kf為聲剛度矩陣；Mf為聲質(zhì)量矩陣。使用有限元法求得聲特征向量p及圓頻率ω，由f=ω/2π即可求得聲振頻率。

1.2 大渦模擬控制方程

在聲渦耦合誘發(fā)固體發(fā)動機不穩(wěn)定燃燒的研究中，漩渦脫落頻率的計算十分關鍵，因而需要采用大渦模擬的方式對發(fā)動機內(nèi)部的湍流狀態(tài)進行高精度求解。大渦模擬的控制方程可通過對原始N-S方程進行空間濾波后得到，以流動參數(shù)f為例，其空間濾波被定義為

(5)

式中：D為積分區(qū)域；g為濾波函數(shù)。對于不可壓流，經(jīng)空間濾波后的N-S方程為

(6)

(7)

2 發(fā)動機壓力振蕩現(xiàn)象分析

2.1 發(fā)動機結構特點

發(fā)動機二脈沖為內(nèi)孔燃燒加局部端面燃燒。二脈沖工作初期發(fā)動機空腔模型如圖1所示。

圖1 發(fā)動機燃燒室空腔模型Fig.1 Numerical model of motor’s combustion chamber

2.2 試驗結果分析

發(fā)動機二脈沖試驗過程壓強-時間曲線如圖2所示。從曲線中可看出發(fā)動機在工作后期出現(xiàn)了較為嚴重的壓力振蕩現(xiàn)象，在1.25 s之前發(fā)動機壓強-時間曲線較為平穩(wěn)，從1.25 s開始出現(xiàn)壓力嚴重的壓力振蕩，直至發(fā)動機工作結束。

將1.25 s后發(fā)動機的振蕩數(shù)據(jù)進行頻域變換后，不同時間段內(nèi)的振蕩主頻如圖3所示，在1～1.25 s時刻內(nèi)燃燒室壓力振蕩比較微弱，其振蕩主頻為320 Hz，對發(fā)動機工作特性幾乎沒有影響；當發(fā)動機工作至1.25 s，振蕩幅值明顯增加，1.25～1.5 s時域內(nèi)的振蕩主頻率為314.7 Hz，隨著發(fā)動機工作時間的增加，壓力振蕩主頻率不斷減小，在發(fā)動機工作末期壓力振蕩主頻率將為284 Hz。

圖2 發(fā)動機試驗過程壓力-時間曲線Fig.2 Pressure versus time for motor in test

圖3 試驗中不同時刻發(fā)動機燃燒室壓力振蕩主頻率Fig.3 Frequencies of motor’s pressure oscillationat various moments in test

2.3 壓力振蕩現(xiàn)象機理分析

圖4 發(fā)動機燃燒室空腔軸向一階聲振型Fig.4 Shape of motor chamber’s first axial acoustic model

根據(jù)藥柱退移形式對不同時刻燃燒室空腔進行建模，利用ANSYS軟件中的acoustic30聲流體單元對不同時刻軸向一階振型的聲頻進行仿真計算，模型表面定義零位移約束。初始時刻燃燒室內(nèi)聲介質(zhì)密度為4.0 kg/m3，平均聲速為1 127 m/s，計算所得初始時刻發(fā)動機軸向一階聲振型如圖4所示，發(fā)動機的軸向一階頻率為351.5 Hz。采用相同方法，計算不同時刻燃燒室空腔的軸向一階振型的聲頻，見表1。從表中可看出發(fā)動機燃燒室軸向一階振型聲頻隨著時間退移呈下降趨勢，與試驗過程中發(fā)動機壓力振蕩主頻趨勢一致，并且計算所得1.25 s后的聲振頻率與試驗結果非常接近，說明發(fā)動機二脈沖工作時出現(xiàn)了以軸向一階聲振頻率為主的不穩(wěn)定燃燒。

表1 不同時刻發(fā)動機軸向一階振型聲頻

利用大渦模擬對不同時刻發(fā)動機燃燒室內(nèi)的流動狀態(tài)進行求解時，需要對網(wǎng)格無關性進行考察。考察結果如圖5所示。當網(wǎng)格數(shù)量達到3.0×105以上時，計算得到的壓力振蕩主頻率已基本保持不變。因此，計算中最終采用的網(wǎng)格量級為3.0×105左右。經(jīng)數(shù)值模擬，發(fā)動機燃燒室內(nèi)的漩渦脫落現(xiàn)象如圖6所示。從圖中可見，二脈沖裝藥的后向臺階是漩渦脫落的主要原因。該漩渦沿軸向不斷向下游發(fā)展，并最終在噴管收斂段與管壁發(fā)生碰撞。同時，由圖4可知，壓力監(jiān)測點1正好處于縱向一階聲振型的波腹，因此其聲壓振幅較大。所以，后續(xù)的數(shù)據(jù)處理主要針對壓力監(jiān)測點1，并將處理得到的燃燒室壓力振蕩主頻率列在表2中。從表2中數(shù)據(jù)可見，發(fā)動機燃燒室內(nèi)漩渦脫落形成的壓力振蕩主頻率同樣隨時間推移不斷下降。

圖5 壓力振蕩主頻率與網(wǎng)格數(shù)量間的關系Fig.5 Relationship between pressure vibration’s frequency and mesh number

圖6 發(fā)動機燃燒室內(nèi)的漩渦脫落Fig.6 Vortex-shedding in motor chamber

時間/s聲頻/Hz時間/s聲頻/Hz0.50338.031.5324.511.00333.672.0295.271.25329.792.5288.49

圖7所示為各個時間段內(nèi)試驗測得的燃燒室壓強振蕩主頻、燃燒室空腔一階軸向聲振頻率和漩渦脫落形成的壓力振蕩主頻三者之間的對比。從圖7中可見，三者十分接近。因此，可以基本確定發(fā)動機工作時發(fā)生的壓力振蕩現(xiàn)象由聲渦耦合引起，即由于燃氣流經(jīng)二脈沖藥柱的后向臺階結構后形成了漩渦脫落，而該漩渦脫落的頻率與發(fā)動機空腔的軸向一階聲振頻率十分接近，從而形成了自激聲振蕩即聲不穩(wěn)定燃燒現(xiàn)象。而且試驗中壓力振蕩現(xiàn)象發(fā)生后，壓力振蕩頻率持續(xù)下降的原因主要是由于空腔聲振頻率和漩渦脫落形成的壓力振蕩主頻同時下降造成的。

圖7 發(fā)動機各組頻率間的比較Fig.7 Comparison for various groups of frequencies

3 擾流環(huán)抑制壓力振蕩現(xiàn)象的機理分析

圖8 不同通徑擾流環(huán)下試驗過程的壓力-時間曲線Fig.8 Pressure versus time for motors using various flow-disturbing rings

為抑制發(fā)動機二脈沖壓力振蕩現(xiàn)象，在發(fā)動機一脈沖與二脈沖之間增加了擾流裝置，并對擾流裝置抑制壓力振蕩現(xiàn)象進行了試驗驗證。試驗結果表明：當擾流環(huán)中心通徑為95 mm時，試驗結果無振蕩；當擾流環(huán)中心通徑為105 mm時，試驗結果基本無振蕩；當擾流環(huán)中心通徑擴大為115 mm時，試驗出現(xiàn)振蕩。不同通徑的擾流環(huán)消振試驗過程的壓力-時間曲線如圖8所示。

針對試驗中采用的3種不同通徑的擾流環(huán)，建立帶擾流環(huán)結構的發(fā)動機燃燒室空腔模型，如圖9所示。對發(fā)動機二脈沖工作時的聲腔頻率和壓力波動頻率進行了計算，分析其抑制壓力振蕩的原理，并探索擾流環(huán)結構自身參數(shù)變化時對壓力振蕩的影響規(guī)律。

圖9 帶擾流環(huán)結構的發(fā)動機燃燒室空腔模型Fig.9 Numerical model of motor’s combustion chamber with flow-disturbing ring

圖10為中心孔徑Φ95 mm擾流環(huán)形成的漩渦結構。由于擾流環(huán)的存在，燃氣在經(jīng)過擾流環(huán)時會形成漩渦，漩渦在往下游移動過程中體積不斷增大。

圖10 中心孔徑Φ95 mm擾流環(huán)形成的漩渦結構Fig.10 Structure of vortex shedding with flow-disturbing ring(Φ95 mm in diameter)

將帶3種尺寸擾流環(huán)結構的發(fā)動機在不同時刻的聲腔頻率、大渦模擬頻率與原發(fā)動機聲腔頻率、大渦模擬結果進行對比，如圖11所示。

圖11 原發(fā)動機與帶不同通徑擾流環(huán) 發(fā)動機各組頻率間的比較Fig.11 Frequency comparison between motors with and without flow-disturbing ring

擾流環(huán)的存在并沒有明顯改變發(fā)動機的軸向一階聲頻，但卻使發(fā)動機工作初期的漩渦脫落頻率得到明顯提高，擾流環(huán)通徑越小發(fā)動機的漩渦脫落頻率越大，即發(fā)動機工作初期的漩渦脫落頻率與發(fā)動機燃燒室空腔的軸向一階聲頻差別越大，從而避免了壓力振蕩現(xiàn)象的出現(xiàn)，這一結論與發(fā)動機的試驗結果相吻合，見表3。

表3 不同擾流環(huán)尺寸下發(fā)動機工作初期(0.5 s時)的頻率對比

4 結束語

本文采用聲腔有限元及大渦模擬技術，對雙脈沖發(fā)動機二脈沖工作過程中的壓力振蕩產(chǎn)生原因及抑制機理進行了研究。研究發(fā)現(xiàn)發(fā)動機二脈沖工作過程中出現(xiàn)的壓力振蕩現(xiàn)象主要由聲渦耦合引起，當燃氣流經(jīng)二脈沖藥柱的后向臺階結構后形成了漩渦脫落，而該漩渦脫落的頻率與發(fā)動機空腔的軸向一階聲振頻率十分接近，從而形成了自激聲振蕩現(xiàn)象；擾流環(huán)對壓力振蕩的抑制作用是由于其提高了發(fā)動機工作初期的漩渦脫落頻率，擾流環(huán)孔徑越小漩渦脫落頻率越高，使該頻率遠離了發(fā)動機空腔的軸向一階聲頻，從而抑制了壓力振蕩的形成。需要說明的是：當聲渦耦合為主要誘因時，文中采用的不穩(wěn)定燃燒分析方法是科學和有效的。所以，它對固體發(fā)動機設計工作具有較好的指導意義，特別是在裝藥設計和不穩(wěn)定燃燒抑制手段的驗證方面。但該方法在使用過程中并沒有考慮兩相流的影響，而凝相顆粒在壓力振蕩的抑制上有時又發(fā)揮著十分重要的作用。因此，后續(xù)研究需要加入這部分內(nèi)容，以使計算模型更加逼近真實情況。