固體火箭發(fā)動機變截面燃燒室縱向一維聲模態(tài)特性研究

史曉鳴， 侯凱宇， 劉陸廣， 高 陽， 夏 鵬， 李海東， 李莎莎， 陸豐瑋

(1. 上海機電工程研究所， 上海 201109； 2. 上海航天技術(shù)研究院， 上海 201109；3. 上海航天動力技術(shù)研究所， 上海 201109； 4. 上海航天精密機械研究所， 上海 201600)

固體火箭發(fā)動機是運載火箭、導(dǎo)彈等細長體飛行器廣泛采用的動力裝置。不穩(wěn)定燃燒是國內(nèi)外固體火箭發(fā)動機研制中屢次遇到的難題[1-3]，其主要表現(xiàn)為發(fā)動機受擾動后，燃燒室內(nèi)周期性壓強振蕩、平均壓強改變，發(fā)動機輸出推力振蕩、平均推力及工作時長改變，并伴有發(fā)動機及飛行器的劇烈振動，嚴重時導(dǎo)致發(fā)動機熄火或超壓爆炸，如圖1所示。對于聲不穩(wěn)定機理的燃燒不穩(wěn)定問題，是固體推進劑燃燒過程與燃燒室內(nèi)聲振過程相互作用的結(jié)果。而固體火箭發(fā)動機往往是細長體構(gòu)型，燃燒室內(nèi)縱向聲模態(tài)頻率遠小于切向和徑向，往往控制著燃燒振蕩的發(fā)展[5]，因此研究縱向聲模態(tài)對于固體火箭發(fā)動機的不穩(wěn)定燃燒具有重要的意義！

經(jīng)典的固體火箭發(fā)動機設(shè)計理論中常將燃燒室簡化為等截面圓柱形，以當?shù)厝細饨橘|(zhì)聲速a和燃燒室長度L，根據(jù)理論解fn=na/2L計算第n階縱向聲模態(tài)固有頻率。但對于翼柱型裝藥、雙脈沖等真實發(fā)動機，在燃燒過程中其內(nèi)部燃燒室截面變化復(fù)雜，理論解往往得不到準確的結(jié)果。文獻[6-8]通過商業(yè)有限元軟件對變截面燃燒室聲模態(tài)的研究表明其縱向固有頻率與等截面理論解存在較大差異，且高階頻率不再是一階基頻的整數(shù)倍。商業(yè)有限元軟件通過對燃燒室三維幾何模型劃分網(wǎng)格進行數(shù)值計算可以精確分析復(fù)雜變截面燃燒室的聲場特性，但每次燃燒室裝藥幾何模型更改設(shè)計后都需要重新劃分網(wǎng)格再進行數(shù)值計算。因此，三維有限元模型不適用于固體火箭發(fā)動機初始方案設(shè)計階段的反復(fù)迭代優(yōu)化及規(guī)律分析，更多的是用于最終設(shè)計狀態(tài)確定后的校核。

圖1 不穩(wěn)定燃燒發(fā)動機燃燒室內(nèi)壓強

考慮到固體火箭發(fā)動機細長體構(gòu)型特征，燃燒室內(nèi)燃氣流動可視作一維加質(zhì)管道內(nèi)流動問題。對一維管道內(nèi)流體聲壓脈動研究是當前熱點。Carmicino等[9]通過求解一維波動方程的解析解，研究了固液混合火箭在點火工作過程中燃燒室軸向聲模態(tài)頻率的變化規(guī)律。Li等[10]推導(dǎo)了環(huán)形燃燒室聲場解析解，研究了非均勻截面和平均流場參數(shù)對聲場特性的影響。Yeddula等[11]提出了變截面管道中的聲波傳播的半解析解。蔡建程等[12]理論分析了不可壓縮流場壓力脈動(偽聲)與聲波方程的區(qū)別，并在CFD(computational fluid dynamics)數(shù)值研究直管與彎管入口處設(shè)置簡諧壓力脈動條件下管內(nèi)非定常流動的特性。趙天泉等[13]設(shè)計了聲學(xué)試驗發(fā)動機，通過試驗研究前封頭結(jié)構(gòu)及潛入空腔對聲穩(wěn)定性的影響規(guī)律。

綜上所述，在固體火箭發(fā)動機設(shè)計過程中，商業(yè)有限元軟件可用于對三維燃燒室聲場特性的精確校核，但在方案階段的反復(fù)迭代優(yōu)化設(shè)計上效率有所不足，且不便于直接分析截面變化對聲模態(tài)影響機理規(guī)律；而一維模型計算效率高，且數(shù)學(xué)模型是個參數(shù)化模型，便于分析及理解其相關(guān)參數(shù)對縱向聲模態(tài)的影響機理及規(guī)律。本文針對固體火箭發(fā)動機燃燒室縱向聲模態(tài)準確高效分析以及變截面對聲模態(tài)影響機理及規(guī)律研究的需求，建立考慮燃燒室截面面積變化的縱向一維聲振方程，有限差分數(shù)值求解縱向聲模態(tài)固有頻率及振型，分析了截面變化對燃燒室縱向聲模態(tài)特性的影響規(guī)律；通過變截面燃燒室模擬試驗件的聲模態(tài)試驗獲取脈動壓力響應(yīng)，基于解析模態(tài)分解和希爾伯特變換辨識縱向聲模態(tài)參數(shù)，驗證了本文數(shù)值方法的準確性；并進一步應(yīng)用于固體火箭發(fā)動機燃燒工作過程中縱向聲模態(tài)的數(shù)值仿真，分析其變化規(guī)律。

1 變截面燃燒室縱向一維聲振方程

固體火箭發(fā)動機燃燒過程中，隨著固體推進劑消耗，燃面不斷退移，燃燒室空間逐漸擴大，且伴有生成燃氣質(zhì)量的注入，其所對應(yīng)的聲振系統(tǒng)是時變的。但聲波在高溫高壓燃氣介質(zhì)內(nèi)的傳播速度遠快于燃面退移和燃氣生成質(zhì)量注入的速度，因此可采用狀態(tài)凝固化的思想，將連續(xù)時變系統(tǒng)在時間上分割為若干個短時時不變系統(tǒng)，提取發(fā)動機工作過程中若干時刻的燃燒室?guī)缀螛?gòu)型及該時刻的燃氣參數(shù)開展聲振分析。

固體火箭發(fā)動機工作過程中某一時刻變截面燃燒室內(nèi)氣體滿足連續(xù)方程

(1)

動量方程

(2)

狀態(tài)方程

(3)

式中：ρ為氣體密度；A=A(x)為燃燒室截面面積；v為聲波脈動速度；p為脈動壓力；a為聲速。

(4)

式(4)對t求導(dǎo)得

(5)

式(2)兩端乘上截面面積A后對x求導(dǎo)得

(6)

(7)

(8)

(9)

2 有限差分法求解縱向一維聲振模態(tài)

對于縱向長度為L的燃燒室，沿縱向均布N個截面，將燃燒室均分為N-1段，每段長Δx=L/(N-1)，第i(i=1～N)個截面上的面積及脈動壓力分別記為Ai，pi。

在第2～第N-1個截面上均滿足式(8)，由此可得N-2個方程

(10)

(11)

(12)

式(12)代入式(10)、式(11)后聯(lián)立得矩陣形式的振動方程

(13)

式中：p={p1,p2,…,pN}T為脈動壓力列向量；M，K分別為聲質(zhì)量、聲剛度矩陣。解廣義特征值(K-λM)·φ=0得聲振模態(tài)。

(14)

(15)

圖2 截面積階梯突變

以及脈動速度與脈動壓強的關(guān)系

(16)

可得

(17)

式(17)表明，脈動壓強對縱向坐標的偏導(dǎo)數(shù)在面積階梯突變截面的兩側(cè)不連續(xù)，截面兩側(cè)的偏導(dǎo)數(shù)與兩側(cè)的面積成反比。

(18)

式(18)代入式(17)得

4pi+1-pi+2)

(19)

式(19)替換式(13)中的第i行即可。

3 試驗驗證及聲模態(tài)辨識方法

真實固體火箭發(fā)動機點火試車試驗成本巨大，且搭載高溫高壓燃氣工況下聲模態(tài)試驗具有危險性。為驗證本文方法的準確性，設(shè)計了若干截面為圓形但截面面積沿縱向變化的燃燒室模擬試驗件，如圖3所示。包括截面半徑線性變化、截面積指數(shù)變化、截面積階梯突變、局部截面收縮等4組工況，用于模擬翼柱型裝藥、雙脈沖等真實固體火箭發(fā)動機中可能存在的各類截面變化情況。試驗件材料為2A12鋁合金，加工過程中，各試驗件一端已一體成型密封，另一端預(yù)留開口用于安裝喇叭，沿外表面一條母線鉆麥克風(fēng)傳感器安裝孔，孔的位置及數(shù)量見圖3，截面半徑線性變化、截面積指數(shù)變化工況均為3個孔，截面積階梯突變、局部截面收縮工況均為5個孔，麥克風(fēng)傳感器編號以最靠近開口端的位置為“1#”，依次遞增類推。

(a) 截面半徑線性變化

(c) 截面積階梯突變

進而開展了一個標準大氣壓下常溫大氣介質(zhì)(聲速a=343 m/s)的縱向聲模態(tài)試驗。試驗系統(tǒng)及現(xiàn)場照片，如圖4所示。模擬試驗件放置于海綿墊上，在開口端安裝HiVi-Q3喇叭并嵌入O型圈密封，沿外表面母線預(yù)留孔安裝Microtech Gefell MM210麥克風(fēng)傳感器并以硅膠密封，Lms.SCR05數(shù)據(jù)采集器輸出猝發(fā)隨機信號，通過DPM200A功率放大器驅(qū)動喇叭對試驗件內(nèi)氣體介質(zhì)施加聲激勵，數(shù)據(jù)采集器同時采集麥克風(fēng)傳感器感受到的脈動壓力信號。

圖4 試驗系統(tǒng)及現(xiàn)場

試驗過程中喇叭實際對試驗件內(nèi)氣體介質(zhì)的聲激勵輸入不可測，只通過麥克風(fēng)傳感器獲取了聲振脈動壓強響應(yīng)。本文采用基于希爾伯特變換的解析模態(tài)分解和振動系統(tǒng)辨識方法從響應(yīng)數(shù)據(jù)中識別聲模態(tài)參數(shù)。

對猝發(fā)隨機激勵下產(chǎn)生的自由衰減脈動壓強信號FFT(fast Fourier transform)后得聲振系統(tǒng)頻譜圖，如圖5所示。圖5中前3階諧振峰清晰易識別，表明本文對模擬試驗件采用猝發(fā)隨機激勵產(chǎn)生自由衰減脈動壓強信號的信噪比良好。以此確定各階聲振諧振頻率ωi,i=(1,2,…,n)。應(yīng)用解析模態(tài)分解方法[14]，選取若干二分頻率ωbi∈(ωi,ωi+1),i=(1,2,…,n-1)，將脈動壓力信號x(t)表示為若干個單頻信號xi(t)的疊加

xi(t)=si(t)-si-1(t)，

s0(t)=0，

si(t)=sin(ωbit)H[x(t)cos(ωbit)]-

cos(ωbit)H[x(t)sin(ωbit)],

i=1,2,…,n-1

(20)

式中：H·為希爾伯特變換；si(t)為頻率小于ωbi∈(ωi,ωi+1),i=(1,2,…,n-1)的信號，即低通信號；xi(t)即為本征模態(tài)信號。典型的自由衰減脈動壓強信號經(jīng)解析模態(tài)分解后的本征模態(tài)信號，如圖6所示。各麥克風(fēng)傳感器位置上的同一階本征模態(tài)信號周期相同，振幅及相位不同，對其進一步分析即可識別模態(tài)參數(shù)。

對每一階本征模態(tài)信號xi(t)應(yīng)用基于希爾伯特變換的振動系統(tǒng)辨識方法[15-16]識別模態(tài)參數(shù)，系統(tǒng)第i階固有圓頻率ωi(t)表達式為

(21)

式中：A(t)為瞬態(tài)幅值或稱包絡(luò)線；ω(t)為信號的瞬態(tài)頻率。

圖5 聲振系統(tǒng)頻譜圖

(a) 1階

(b) 2階

(c) 3階

對于聲阻尼較小的系統(tǒng)，系統(tǒng)第i階聲模態(tài)在第a個傳感器測點上的振型可以表示為

Sia=Aai(t)cosψi/Abi(t)

(22)

式中：Aai(t)，Abi(t)分別為a，b兩個測點上的第i階本征模態(tài)信號的包絡(luò)線；ψi為a，b兩個測點上的第i階本征模態(tài)信號的相位差。

4 結(jié)果與討論

4.1 變截面燃燒室模擬試驗件縱向聲模態(tài)

變截面燃燒室模擬試驗件縱向聲模態(tài)頻率采用本文方法計算、有限元軟件ABAQUS三維實體模型仿真以及試驗結(jié)果及其對比誤差，如表1、表2所示。從表1、表2中可見：本文方法與有限元仿真及試驗的結(jié)果比較均非常接近，截面半徑線性變化、截面積指數(shù)變化兩種工況下的計算與試驗結(jié)果誤差略微大于截面積階梯突變、局部截面收縮兩種工況下的計算與試驗結(jié)果誤差。這是由于麥克風(fēng)傳感器安裝時，其咪頭需略突出于內(nèi)表面以充分感受到壓強脈動，從而引起該安裝截面上的截面積局部小幅度突變，導(dǎo)致試驗與仿真結(jié)果之間存在誤差，如圖7所示；但在截面積階梯突變、局部截面收縮兩種工況下，試驗件本身存在較大幅度的局部截面積突變，故麥克風(fēng)安裝引起的局部截面積變化對試驗結(jié)果的影響反而不明顯。

從表1中可見，變截面試驗件第n階(n>1)縱向聲模態(tài)頻率fn與 1階基頻f1之間的關(guān)系fnf1≠n，而是n-1

前3階縱向聲模態(tài)振型S1，S2，S3如圖8所示。圖8中左端對應(yīng)封閉端，右端對應(yīng)開口安裝喇叭端。本文約定以1#麥克風(fēng)傳感器所在位置的振型為1進行歸一化，該位置在圖8各個圖中均為最右側(cè)的試驗測點。變截面工況下，縱向聲模態(tài)的奇數(shù)階振型關(guān)于縱向長度的中點不再保持反對稱，偶數(shù)階振型關(guān)于縱向長度的中點不再保持對稱：以與燃燒振蕩密切相關(guān)的一階聲模態(tài)為例，截面半徑線性漸變、截面積指數(shù)變化、截面積階梯突變等變截面工況的縱向聲模態(tài)振型的最大幅值均出現(xiàn)在截面積較小的一側(cè)端面，同時其振型節(jié)點位置與等截面振型節(jié)點相比，不再處于縱向長度的中點位置，而是從中點向截面積較大的一側(cè)移動；局部截面收縮變截面工況的縱向聲模態(tài)振型的最大幅值出現(xiàn)在中間收縮段的一側(cè)端面，同時其振型節(jié)點位置與等截面振型節(jié)點相比，不再處于中點位置，而是從中點向局部收縮段的一側(cè)移動。這是由于在變截面模擬試驗件內(nèi)為同一氣體介質(zhì)，氣體介質(zhì)體積模量K是均勻分布的，但截面積A(x)沿著縱向變化，故各個截面上的氣體介質(zhì)拉壓剛度K·A(x)沿著縱向分布是不同的，截面積較大區(qū)域其拉壓剛度較大，脈動聲壓的振型幅值會較小，由此必然造成節(jié)點兩側(cè)的振型不再保持對稱或者反對稱。

表1 縱向聲模態(tài)頻率

表2 縱向聲模態(tài)頻率結(jié)果誤差

圖7 麥克風(fēng)傳感器咪頭安裝示意圖

4.2 固體火箭發(fā)動機燃燒工作過程中縱向聲模態(tài)

進一步將本文方法應(yīng)用于某內(nèi)孔型裝藥真實固體火箭發(fā)動機燃燒工作過程中的縱向聲模態(tài)分析。取整個工作過程中歸一化時間tn=0，1/3，1/2，2/3，1共5個時刻的燃燒室?guī)缀螛?gòu)型開展縱向聲模態(tài)分析。前3階縱向聲頻率變化規(guī)律，如圖9所示。在發(fā)動機工作的初始階段，其縱向聲模態(tài)頻率明顯小于等截面狀態(tài)下的頻率；隨著發(fā)動機燃燒工作直至熄火，縱向聲模態(tài)頻率逐漸單調(diào)上升趨向于等截面狀態(tài)下的頻率，直至最終接近等截面狀態(tài)下的頻率。

圖9 固體火箭發(fā)動機燃燒工作過程中縱向聲模態(tài)頻率

圖10 固體火箭發(fā)動機燃燒工作過程中燃燒室截面積變化

5 結(jié) 論

(1) 本文推導(dǎo)了考慮燃燒室截面面積變化的縱向一維聲振方程，有限差分數(shù)值求解縱向聲模態(tài)固有頻率及振型；設(shè)計了變截面燃燒室模擬試驗件，以聲模態(tài)試驗驗證了本文數(shù)值方法的準確性。

(2) 變截面燃燒室模擬件試驗表明：變截面工況下，縱向聲模態(tài)頻率第n階(n>1)縱向聲模態(tài)頻率fn與 1階基頻f1之間的關(guān)系fnf1≠n，而是n-1

(3) 對于內(nèi)孔翼柱型裝藥結(jié)構(gòu)的固體火箭發(fā)動機，在發(fā)動機初始工作階段，截面變化最為劇烈，工程中常用的將燃燒室簡化為等截面圓柱分析方法誤差較大，此時必須考慮截面變化的影響，從而得到更為準確的聲模態(tài)。