過載條件下固體發(fā)動機聲能共振規(guī)律研究①

宋儒儒,張翔宇,甘曉松,趙天泉

(1.中國航天科技集團有限公司四院四十一所,西安 710025;2.中國航天科技集團有限公司第四研究院,西安 710025)

0 引言

固體火箭發(fā)動機用于大型運載火箭助推器和各類導彈武器已超過半個世紀,但不穩(wěn)定燃燒現(xiàn)象在各種尺度的戰(zhàn)術(shù)導彈發(fā)動機及分段式發(fā)動機上時常發(fā)生,是發(fā)動機研制中不可忽略的問題[1-6]。發(fā)動機的不穩(wěn)定燃燒本質(zhì)上是其內(nèi)部各種增益/阻尼因素之間相互作用導致的壓強周期性振蕩。粒子阻尼作為其中的一項阻尼因素,占總阻尼的1/3以上,遠遠大于結(jié)構(gòu)阻尼和氣相阻尼[1]。

粒子阻尼理論研究中使用最廣的是微粒松弛理論,當惰性微粒的動力弛豫時間與聲振蕩的特征時間相等時,惰性微粒的阻尼效應最大。CAI等[7]基于歐拉-拉格朗日法對固體火箭發(fā)動機內(nèi)兩相流與聲能的交換機理進行了數(shù)值模擬,充分考慮了聲場、湍流散布、粒子碰撞與聚合等復雜流場條件下氣相和粒子相之間的耦合作用,指出粒子的動力馳豫時間、熱馳豫時間和聲場的特征時間之間的關(guān)系對發(fā)動機內(nèi)流場兩相流的耦合起著非常重要的作用。西北工業(yè)大學劉洋、何國強等[8]提出了一種可收集燃燒室中凝聚態(tài)粒子的方法,設計了粒子收集試驗裝置,通過改變試驗裝置的收斂半角和試驗狀態(tài)參數(shù)可以模擬真實發(fā)動機中的粒子聚集狀態(tài)。彭小波等[9]改進了燃燒室中凝相燃燒產(chǎn)物的試驗收集裝置,并開展了不同聚集狀態(tài)下凝相燃燒產(chǎn)物的燃燒特性分析。金秉寧等[10]針對不同初始鋁粉粒度的含鋁復合推進劑,對其燃燒產(chǎn)物進行了粒子阻尼特性和鋁粒子分布燃燒響應特性的試驗研究,得到了關(guān)于粒子阻尼的大小及預估、粒子對分布燃燒的影響方面的規(guī)律。

隨著對粒子阻尼的深入研究,可以從發(fā)動機飛行狀態(tài)變化和燃燒室聲腔聲壓振蕩等動態(tài)的角度理解粒子阻尼。固體火箭在做機動飛行時,橫向過載往往會改變?nèi)紵覂?nèi)凝相粒子的空間分布[11-12],從而影響粒子的阻尼特性,甚至會導致不穩(wěn)定燃燒。西北工業(yè)大學劉佩進課題組[13]對此展開研究,研究結(jié)果表明在大過載條件下,燃燒室內(nèi)大尺寸凝相顆粒會在極短的時間內(nèi)在局部形成聚集帶和真空帶,導致粒子阻尼減小,這對發(fā)動機穩(wěn)定性是十分不利的。同時,在極限環(huán)振蕩的發(fā)展過程中,由于聲對粒子的操縱作用,可能導致粒子沿發(fā)動機縱向方向重新分布,引起粒子阻尼特性的改變,并參與壓強振蕩非線性增長過程中的增益與阻尼博弈[14]。北京理工大學王寧飛團隊[15-16]使用脈沖衰減法對粒子阻尼進行數(shù)值仿真研究,結(jié)果表明,由脈沖衰減法計算所得的結(jié)果與T型燃燒器的實驗結(jié)果基本吻合,驗證了該仿真方法在計算粒子阻尼時的合理性與有效性。

本文從聲學角度出發(fā),用聲場介質(zhì)的聲衰減系數(shù)代替粒子阻尼,探究固體火箭發(fā)動機在不同飛行過載條件下聲學特性的變化規(guī)律,以聲學響應傳遞函數(shù)來評估發(fā)動機的聲不穩(wěn)定性。

1 粒子阻尼理論

粒子阻尼作用是由弛豫現(xiàn)象引起的。由于氣體的粘性和粒子的不可壓縮性,使粒子與氣體達到相同的速度、溫度需要一定的時間,所以存在動力弛豫過程與熱力弛豫過程。

定義動力弛豫時間τv,即粒子速度達到與氣體速度平衡所需要的時間。

(1)

式中ms為單個粒子的質(zhì)量,kg;μ為氣體的動力粘度系數(shù),kg/(m·s);σ為粒子半徑,μm;ρs為粒子自身密度,kg/m3。

定義熱力弛豫時間τt,是粒子溫度達到與氣體溫度平衡所需要的時間。

(2)

式中Pr為普朗特數(shù),Pr=μcp/λg;cs為粒子的比熱容,J/(kg·K);cp為氣體比定壓熱容,J/(kg·K)。

用壓強振幅的時間衰減率(即衰減系數(shù))αp來表示粒子阻尼的大小,粒子阻尼系數(shù)可表示為

(3)

以上各式中動力粘度系數(shù)μ由Sutherland定律確定:

式中μ0為參考粘度,1.716×10-5kg/(m·s);T為燃氣靜溫;T0為參考溫度,273.11 K;S為Sutherland常數(shù),110.56 K。

2 聲能共振數(shù)值計算方法

2.1 總體思路

本文通過將粒子阻尼轉(zhuǎn)化為聲場介質(zhì)的聲衰減系數(shù),以聲學響應傳遞函數(shù)評估發(fā)動機聲不穩(wěn)定性,研究過程中,不考慮其他阻尼因素影響,思路如下:

(1)對發(fā)動機幾何模型劃分計算域,并進行聲腔模態(tài)提取及噴管阻尼計算,以噴管阻尼計算相應粒子阻尼范圍,確定凝相粒子粒徑與濃度的初始參數(shù);

(2)進行兩相流計算,并統(tǒng)計各個計算域的粒子濃度,將各計算域濃度結(jié)果帶入式(3)計算相應粒子阻尼值;

(3)用聲場介質(zhì)的聲衰減系數(shù)代替粒子阻尼值,表征粒子阻尼對聲腔聲學特性的影響;

(4)在各計算域添加相應聲衰減系數(shù),計算發(fā)動機聲腔的聲學特性,分析聲學響應傳遞函數(shù)的變化規(guī)律。

2.2 計算模型

本文選用某發(fā)動機工作后期的簡化模型作為研究對象。由于噴管喉部下游為超聲速流動,聲波不再反射,因此在進行仿真分析時選取噴管喉部上游部分結(jié)構(gòu),幾何模型及主要參數(shù)如圖1所示。為統(tǒng)計過載條件下發(fā)動機內(nèi)凝相粒子濃度分布,對模型進行計算域劃分,并參照笛卡爾坐標系方向?qū)Ω鲄^(qū)域編號,如圖2所示。為便于表述,將所有計算域劃分為4部分:第一部分為坐標系(+Y+Z)區(qū)域,包含1～37號計算域,記作A;第二部分為坐標系(+Y-Z)區(qū)域,包含2～38號計算域,記作B;第三部分為坐標系(-Y-Z)區(qū)域,包含3～39號計算域,記作C;第四部分為坐標系(-Y+Z)區(qū)域,包含4～40號計算域,記作D。

圖1 某發(fā)動機工作后期幾何構(gòu)型

圖2 不同計算域編號

提取聲腔模態(tài),模態(tài)分布如圖3所示,模態(tài)數(shù)值見表1。參考穩(wěn)態(tài)波衰減法計算該模型的噴管阻尼,在發(fā)動機頭部添加單極子聲源,聲源頻率設置為與聲腔固有頻率相同,本文取162.094 Hz,幅值為0.1 kg/s,并在聲腔內(nèi)建立起穩(wěn)定的駐波后關(guān)閉聲源。聲腔內(nèi)壓力振蕩幅值以指數(shù)形式衰減,將衰減過程的振蕩幅值繪制在半對數(shù)坐標系中取包絡線并進行線性擬合,所得擬合直線斜率即為衰減系數(shù),聲腔內(nèi)聲壓衰減過程以及擬合直線斜率如圖4所示。

(a)Sound pressure attenuation process

表1 聲腔各階模態(tài)值

(a)First natural frequency

對聲壓衰減過程進行直線擬合,斜率k=-25.22,即噴管阻尼αN=25.22 s-1,分析發(fā)動機內(nèi)部各種聲能阻尼占比[1],αN∶αp≈3∶2,因此該發(fā)動機模型中粒子阻尼在16 s-1左右。經(jīng)反復校核,在粒徑為48 μm,濃度為15%時,通過公式(3)計算所得的粒子阻尼為16.8 s-1,符合兩種阻尼在總阻尼中的占比,可進行后續(xù)計算。

2.3 計算工況參數(shù)

氣相控制方程采用N-S方程,顆粒相運動軌跡則由拉格朗日方程確定,獲得粒子速度和濃度的空間分布。穩(wěn)態(tài)氣相流場計算采用質(zhì)量流量入口,入口流量為6.072 kg/s,介質(zhì)溫度為3235.88 K,出口為壓力出口邊界,喉部出口靜壓為4 035 328 Pa,靜溫為300 K,壁面選擇無滑移邊界條件,使用基于壓力的求解器,湍流模型選用SSTk-ε模型,對控制方程進行二階迎風格式離散,計算結(jié)果如圖5所示。

(a)Gas phase pressure distribution

由固體發(fā)動機原理[17]可知,燃燒室末端燃氣速度系數(shù)一般為0.2～0.5,即燃氣速度為200～500 m/s。以長度為3 m的發(fā)動機為例,據(jù)此估算粒子在燃燒室內(nèi)最長停留時間為30 ms,因此本文計算時間取30 ms。凝相粒子的初始速度取為氣相入口速度的0.4倍[18],壁面條件為彈性反彈,凝相粒子粒徑、濃度及過載參數(shù)見表2。

表2 凝相粒子初始參數(shù)

2.4 粒子阻尼的轉(zhuǎn)換

聲場中介質(zhì)的阻尼效應通常使用聲衰減系數(shù)表征,即在聲速表達式中賦加虛部值,如式(4)所示。

(4)

通過擬合聲衰減系數(shù)Y與粒子阻尼系數(shù)αp的關(guān)系,可以將粒子阻尼轉(zhuǎn)化為對應的聲衰減系數(shù)。對Y分別取10組值,利用穩(wěn)態(tài)波衰減法計算聲腔內(nèi)阻尼值,計算方法同2.2節(jié),計算所得阻尼值減去噴管阻尼值即為粒子阻尼αp,擬合結(jié)果如圖6所示。

圖6 Y與αp 的關(guān)系擬合

粒子阻尼αp與聲衰減系數(shù)Y關(guān)系式如式(5):

Y=0.152 6αp+0.253 128

(5)

2.5 聲學仿真計算

采用直接頻率響應分析方法,在發(fā)動機頭部添加球狀聲源,聲源幅值為2000 Pa,并在每個計算域中心位置添加一個聲壓頻率監(jiān)測點,同時定義聲學響應傳遞函數(shù)Pf來表征壓力振蕩程度,見式(6),依據(jù)該參數(shù)來評估發(fā)動機出現(xiàn)聲不穩(wěn)定的趨勢。

(6)

式中p為測點聲壓值,Pa;ps為聲源聲壓值,Pa。

3 結(jié)果及分析

3.1 各計算域粒子濃度分布規(guī)律

設置軸向過載為15g,橫向過載分別為3g、6g、9g、12g、15g,過載施加時間同2.3節(jié)仿真計算時間一致,即從0 s開始直至仿真計算結(jié)束。探究30 ms時刻下粒子空間分布規(guī)律以及各工況下聲學響應情況,同時以未加過載即橫向過載與軸向過載均為0g的工況作為對照,各工況下粒子濃度的空間分布如圖7、圖8所示。

(a)Transverse overload=0g (b)Transverse overload=3g

(a)Transverse overload=3g (b)Transverse overload=6g (c)Transverse overload=9g

采用側(cè)向加質(zhì),軸向及橫向過載方向分別與X、Y軸正方向保持一致。結(jié)合圖7～圖9,在施加軸向過載的工況下發(fā)動機內(nèi)所有粒子的速度方向與X軸相同,不斷向噴管喉部運動。發(fā)動機頭部基本不存在粒子軌跡,形成粒子真空區(qū),噴管喉部的粒子不斷聚集,形成粒子聚集區(qū)。

在橫向過載的影響下,粒子不斷沿Y軸矢量方向在發(fā)動機近壁面區(qū)域聚集,如圖8所示。并隨著橫向過載由3g不斷增加到15g的過程中,粒子聚集區(qū)域中心沿Y軸的偏移量逐漸增大。當橫向過載加至12g和15g時,粒子聚集區(qū)域縮小,粒子空間分布更加集中,且聚集區(qū)中心沿Y軸偏移較大,同時在Y軸負方向發(fā)動機近壁面區(qū)域形成較大的粒子真空區(qū),如圖8(d)、(e)所示。

在橫向、軸向過載的共同作用下,發(fā)動機噴管附近A、B區(qū)域粒子濃度明顯高于C、D區(qū)域,如圖9所示。

(a)Domain 1～37 (b)Domain 2～38

3.2 網(wǎng)格無關(guān)性驗證

為檢驗仿真計算的合理性,對模型繪制800 000、1 000 000、1 200 000數(shù)目的網(wǎng)格,進行網(wǎng)格無關(guān)性驗證,網(wǎng)格模型如圖10所示。以表2中case 6為例,分別計算各套網(wǎng)格模型中Domain 1～37計算域的粒子阻尼,結(jié)果如圖11所示。

(a)80×104

由圖11可知,不同網(wǎng)格數(shù)目下,粒子阻尼值隨計算域監(jiān)測點的變化趨勢基本一致,即阻尼計算結(jié)果對網(wǎng)格依賴性極小。在計算誤差允許范圍內(nèi),可證明網(wǎng)格數(shù)量對數(shù)值計算結(jié)果無影響。

圖11 不同網(wǎng)格數(shù)目下的阻尼值

3.3 粒子濃度分布對聲學響應的影響

將3.1節(jié)各區(qū)域計算所得的粒子濃度帶入式(3)計算對應粒子阻尼值,計算所得阻尼值結(jié)果帶入式(5),可計算相應聲衰減系數(shù),將其添加至聲場中進行聲學響應計算。

對未加過載的工況進行頻率響應分析,取37、38、39、40等四個區(qū)域中心添加監(jiān)測點,繪制頻率響應曲線如圖12所示?？梢钥闯鏊膫€監(jiān)測點的頻率響應曲線重合,計算結(jié)果一致,且發(fā)動機內(nèi)模態(tài)頻率分布以一階頻率為主,故后續(xù)對不同工況下各監(jiān)測點聲學響應的分析主要在一階頻率下進行。

圖12 零過載工況下各測點頻率響應曲線

發(fā)動機模型具有很好的對稱性,故只分析A、D區(qū)域測點的聲學響應。不同工況下各計算域一階振蕩振幅分布如圖13所示,發(fā)動機頭部和噴管喉部測點的聲壓振幅最大,發(fā)動機中間位置振幅較小,符合一階軸向振型分布特點,且發(fā)動機頭部聲學響應略高于噴管喉部。

(a)Domain 1～37 (b)Domain 4～40

對發(fā)動機頭部與噴管喉部區(qū)域進行聲學響應分析,如圖14所示,可見所有監(jiān)測點的聲學響應隨橫向過載的增加不斷增大。圖14(a)是Domain 1和Domain 37的監(jiān)測點聲學響應分布情況,Domain 37為粒子聚集區(qū),并隨著橫向過載的不斷增加,該區(qū)域的粒子濃度不斷增加,但聲學響應也呈不斷增大的趨勢,并未出現(xiàn)由于粒子聚集使該區(qū)域的聲學響應減小,這是因為橫向過載導致粒子的整體空間分布極其不均勻,對聲振蕩的阻尼作用減小,導致聲學響應增大。對比Domain 1和Domain 37的聲學響應,橫向過載相同時,Domain 1的聲學響應高于Domain 37,符合粒子阻尼理論,即粒子濃度增加,粒子阻尼值增加,對聲不穩(wěn)定的抑制作用越強。圖14(b)所示為Domain 4和Domain 40的監(jiān)測點聲學響應分布情況,在橫向過載由3g不斷增加到15g的過程中,聲學響應呈現(xiàn)出不斷增加的趨勢,且相同過載情況下Domain 4的聲學響應高于Domain 40。

(a)Domain 1 and Domain 37 (b)Domain 4 and Domain 40

4 結(jié)論

(1)在粒子阻尼不斷減小時,聲腔內(nèi)的聲學響應傳遞函數(shù)不斷增大,發(fā)動機出現(xiàn)不穩(wěn)定燃燒現(xiàn)象的概率增加。

(2)在凝相粒子粒徑保持一致且計算頻率為一階固有頻率的前提下,隨著橫向過載的不斷增加,發(fā)動機內(nèi)粒子沿過載方向明顯聚集,并且發(fā)動機內(nèi)各測點處的聲學響應值隨橫向過載的增大呈增大趨勢,表明大橫向過載工況下,發(fā)動機內(nèi)粒子的阻尼作用減小,這對發(fā)動機燃燒穩(wěn)定性是極其不利的。

(3)同一過載工況下,發(fā)動機頭部跟噴管喉部的聲學響應最大,中間位置聲學響應最小,符合一階壓力振蕩分布特點。同時頭部聲學響應值高于噴管喉部,這是因為軸向過載使發(fā)動機內(nèi)粒子向噴管附近聚集,增大了對聲振蕩的抑制作用。