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

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

0 引言

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

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

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

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

1 粒子阻尼理論

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

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

(1)

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

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

(2)

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

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

(3)

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

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

2 聲能共振數值計算方法

2.1 總體思路

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

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

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

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

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

2.2 計算模型

本文選用某發(fā)動機工作后期的簡化模型作為研究對象。由于噴管喉部下游為超聲速流動,聲波不再反射,因此在進行仿真分析時選取噴管喉部上游部分結構,幾何模型及主要參數如圖1所示。為統(tǒng)計過載條件下發(fā)動機內凝相粒子濃度分布,對模型進行計算域劃分,并參照笛卡爾坐標系方向對各區(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ā)動機工作后期幾何構型

圖2 不同計算域編號

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

(a)Sound pressure attenuation process

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

(a)First natural frequency

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

2.3 計算工況參數

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

(a)Gas phase pressure distribution

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

表2 凝相粒子初始參數

2.4 粒子阻尼的轉換

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

(4)

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

圖6 Y與αp 的關系擬合

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

Y=0.152 6αp+0.253 128

(5)

2.5 聲學仿真計算

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

(6)

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

3 結果及分析

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

設置軸向過載為15g,橫向過載分別為3g、6g、9g、12g、15g,過載施加時間同2.3節(jié)仿真計算時間一致,即從0 s開始直至仿真計算結束。探究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

采用側向加質,軸向及橫向過載方向分別與X、Y軸正方向保持一致。結合圖7～圖9,在施加軸向過載的工況下發(fā)動機內所有粒子的速度方向與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 網格無關性驗證

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

(a)80×104

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

圖11 不同網格數目下的阻尼值

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

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

對未加過載的工況進行頻率響應分析,取37、38、39、40等四個區(qū)域中心添加監(jiān)測點,繪制頻率響應曲線如圖12所示?？梢钥闯鏊膫€監(jiān)測點的頻率響應曲線重合,計算結果一致,且發(fā)動機內模態(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ū)域的粒子濃度不斷增加,但聲學響應也呈不斷增大的趨勢,并未出現由于粒子聚集使該區(qū)域的聲學響應減小,這是因為橫向過載導致粒子的整體空間分布極其不均勻,對聲振蕩的阻尼作用減小,導致聲學響應增大。對比Domain 1和Domain 37的聲學響應,橫向過載相同時,Domain 1的聲學響應高于Domain 37,符合粒子阻尼理論,即粒子濃度增加,粒子阻尼值增加,對聲不穩(wěn)定的抑制作用越強。圖14(b)所示為Domain 4和Domain 40的監(jiān)測點聲學響應分布情況,在橫向過載由3g不斷增加到15g的過程中,聲學響應呈現出不斷增加的趨勢,且相同過載情況下Domain 4的聲學響應高于Domain 40。

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

4 結論

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

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

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