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

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

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

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

圖1 不穩(wěn)定燃燒發(fā)動(dòng)機(jī)燃燒室內(nèi)壓強(qiáng)

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

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

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

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

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

(1)

動(dòng)量方程

(2)

狀態(tài)方程

(3)

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

(4)

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

(5)

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

(6)

(7)

(8)

(9)

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

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

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

(10)

(11)

(12)

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

(13)

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

(14)

(15)

圖2 截面積階梯突變

以及脈動(dòng)速度與脈動(dòng)壓強(qiáng)的關(guān)系

(16)

可得

(17)

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

(18)

式(18)代入式(17)得

4pi+1-pi+2)

(19)

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

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

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

(a) 截面半徑線性變化

(c) 截面積階梯突變

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

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

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

對(duì)猝發(fā)隨機(jī)激勵(lì)下產(chǎn)生的自由衰減脈動(dòng)壓強(qiáng)信號(hào)FFT(fast Fourier transform)后得聲振系統(tǒng)頻譜圖，如圖5所示。圖5中前3階諧振峰清晰易識(shí)別，表明本文對(duì)模擬試驗(yàn)件采用猝發(fā)隨機(jī)激勵(lì)產(chǎn)生自由衰減脈動(dòng)壓強(qiáng)信號(hào)的信噪比良好。以此確定各階聲振諧振頻率ωi,i=(1,2,…,n)。應(yīng)用解析模態(tài)分解方法[14]，選取若干二分頻率ωbi∈(ωi,ωi+1),i=(1,2,…,n-1)，將脈動(dòng)壓力信號(hào)x(t)表示為若干個(gè)單頻信號(hào)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)的信號(hào)，即低通信號(hào)；xi(t)即為本征模態(tài)信號(hào)。典型的自由衰減脈動(dòng)壓強(qiáng)信號(hào)經(jīng)解析模態(tài)分解后的本征模態(tài)信號(hào)，如圖6所示。各麥克風(fēng)傳感器位置上的同一階本征模態(tài)信號(hào)周期相同，振幅及相位不同，對(duì)其進(jìn)一步分析即可識(shí)別模態(tài)參數(shù)。

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

(21)

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

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

(a) 1階

(b) 2階

(c) 3階

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

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

(22)

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

4 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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