金 烜，沈赤兵，林 森，楊林濤

(1. 國防科技大學 空天科學學院 高超聲速沖壓發(fā)動機技術重點實驗室, 湖南 長沙 410073; 2. 西安衛(wèi)星測控中心, 陜西 西安 710043)

高頻燃燒不穩(wěn)定是液體火箭發(fā)動機研制過程中極具挑戰(zhàn)性的課題[1]。發(fā)生高頻燃燒不穩(wěn)定時，燃燒室壁面熱流增大導致發(fā)動機在幾分之一秒內被燒毀，高幅值的壓力振蕩甚至會引起發(fā)動機瞬間爆炸。目前，關于高頻燃燒不穩(wěn)定性的維持機理尚未達成共識，預測燃燒穩(wěn)定性依然是新型發(fā)動機研制過程中最困難的一環(huán)[2]，一個原因就是噴注器采用不同的推進劑組合，引起霧化、蒸發(fā)、混合與燃燒等過程發(fā)生變化。工程上往往采用大量全尺寸試車來克服高頻燃燒不穩(wěn)定，但試錯法的燃燒試驗成本和迭代設計成本高昂，因此能在發(fā)動機預研階段采用理論分析、數(shù)值仿真和試驗驗證等多種手段來探索不穩(wěn)定機理有著巨大的經(jīng)濟價值和時間效益。

國內外對液體火箭發(fā)動機非線性不穩(wěn)定燃燒機理[3-5]和各類抑制裝置[6-8]已開展了系統(tǒng)的研究。同時隨著推進劑無毒化逐漸成為未來的主流，常規(guī)的同軸式噴注器和撞擊式噴注器分別采用液氧/液氫等低溫推進劑組合與液氧/煤油等常溫推進劑組合，已開展大量燃燒穩(wěn)定性研究并廣泛應用于新型運載火箭[3, 9-10]。

起源于20世紀50年代的針栓式噴注器具有結構簡單、面積可調和燃燒穩(wěn)定等諸多優(yōu)點[11]，若采用液氧/甲烷雙組元推進劑，可在未來星際探索中實現(xiàn)飛行器就地加注和重復使用[12]。針栓式發(fā)動機燃燒特性關于幾何尺寸和噴注條件的參數(shù)化研究已見諸報道[13-17]，但其中關于其聲學振蕩特性的研究還較少，且局限于燃燒室壓力采集及其頻譜分析、非穩(wěn)態(tài)噴霧燃燒的光學觀測及其模態(tài)分解等分析手段。

相比于試驗和非線性理論分析，數(shù)值模擬振蕩燃燒可獲取更詳細的流場信息。前期的仿真工作[15]說明LOX/GCH4針栓式發(fā)動機工作時燃燒穩(wěn)定性較好，燃燒室壓力以低頻小幅振蕩為主。本文在此基礎上借鑒文獻[9-10]中旋轉齒輪聲學擾動裝置的思路，通過對燃燒流場施加速度擾動來產生橫向聲學振蕩響應，以期在熱試之前預先了解聲學激勵頻率對非穩(wěn)態(tài)噴霧燃燒的影響。

1 數(shù)值物理模型

1.1 物理模型

LOX/GCH4針栓式火箭發(fā)動機樣機采用設計流量為356 g/s的液氧中心式氣液針栓噴嘴(液氧和氣甲烷的設計流量分別為257 g/s和99 g/s)。如圖1(a)和圖1(b)所示，液氧通過針栓中心流道后自針栓頭端的一系列孔呈放射狀徑向噴出，氣甲烷從噴注器集氣腔流經(jīng)針栓外側的環(huán)縫軸向噴出形成氣膜，徑向射流與軸向氣膜呈90°交叉撞擊后推進劑發(fā)生霧化、混合與燃燒，局部動量比為1，相關結構參數(shù)如下：Ls=12 mm,Dp=12 mm,h=0.66 mm,d=0.84 mm,N=12。數(shù)值模擬采用的計算域見圖1(c)，全文坐標原點取燃燒室入口橫截面的中心點，流動方向為x，寬度方向為z，高度方向為y；本文通過矩形燃燒室方案實現(xiàn)用橫向振蕩模擬圓筒形燃燒室內可能出現(xiàn)的復雜切向振蕩，同時矩形燃燒室方案有利于在壁面局部位置開設觀察窗，便于熱試過程中同步觀測噴霧分布和火焰形態(tài)(開窗觀測手段已在針栓式液體火箭發(fā)動機中多次應用[5, 13-14])；矩形燃燒室內部空間寬45 mm，高141 mm，長348 mm，噴管喉部處高收縮至8.2 mm，噴管出口處高擴展至11.9 mm；數(shù)值計算時還在z=0平面上設置多個壓力監(jiān)測點，其中P4和P7位于燃燒室軸線，其余位于壁面，各點對應的x軸坐標如下：-28 mm(P01/P02), 0(P1/P2), 174 mm(P3/P4/P5), 174 mm(P6/P7/P8)。

(a) 針栓式發(fā)動機結構示意圖(a) Schematic of pintle engine

1.2 計算方法和網(wǎng)格劃分

LOX/GCH4針栓式火箭發(fā)動機的兩相噴霧燃燒流動采用Euler-Lagrange方法開展計算，控制方程包含Euler坐標系下的雷諾平均納維-斯托克斯(Reynolds Averaged Navier-Stokes,RANS)方程(氣相)和Lagrange坐標系下的離散顆粒軌道方程(液滴)，兩相之間的耦合通過氣液相互作用源項實現(xiàn)。采用有限體積法離散控制方程，而后采用密度基隱式算法進行求解?？紤]到計算資源有限，計算中采用Jones-Lindstedt四步反應機理結合渦耗散概念模型對GOX/GCH4化學反應過程進行模擬[18]。

計算域構型(見圖1(c))中GCH4的入口邊界給定質量流量和溫度(300 K)；液氧噴注采用離散相模型(Discrete Phase Model, DPM)，以針栓頭端的液氧噴孔面(見圖1(b))作為入口噴射源，顆粒類型為Droplet，所有液滴的噴注方向垂直于噴孔面，并設置液滴的質量流量、入口溫度(80 K)、入口速度(33.48 m/s，對應噴注壓降為1 MPa)和液滴粒徑(10 μm)；出口邊界由數(shù)值計算外推給定。所有壁面邊界采用絕熱非滑移壁面條件，液滴與壁面發(fā)生碰撞后發(fā)生球形反彈。特別說明的是，液滴粒徑的選擇是基于前期的仿真研究[15]，燃燒室壓力振蕩幅度隨著液滴粒徑增大而上升，當液滴粒徑大于60 μm時出現(xiàn)低頻不穩(wěn)定燃燒，進一步增至110 μm以上則發(fā)生熄火。為避免低頻不穩(wěn)定燃燒對聲學振蕩的影響和發(fā)動機性能下降，此處液滴粒徑取10 μm。

由于計算域構型結構不規(guī)則，噴注器集氣腔和環(huán)縫部分采用結構網(wǎng)格，燃燒室和噴管部分采用非結構網(wǎng)格，兩部分在環(huán)縫出口面進行拼接。其中，非結構網(wǎng)格部分采用自下而上的方法生成：先在不同邊界處生成密度不同的面網(wǎng)格，在此基礎上于針栓頭端表面生成的邊界層網(wǎng)格，最后在其余空間里充填局部加密的體網(wǎng)格，從而保證混合反應區(qū)和尾噴管有相對較密的網(wǎng)格。整個計算域包含74萬個網(wǎng)格單元，滿足標準k-ε湍流模型和壁面函數(shù)的使用要求，時間步長取5×10-6s。

關于計算方法及其驗證和網(wǎng)格無關性的詳細介紹參見文獻[15]，同時文獻[5]通過對比封閉方腔內的一維聲學振蕩理論解和二維仿真結果，證明了本文的計算方法適合于聲學振蕩數(shù)值模擬。

1.3 聲學激勵方法

本文的計算順序是先按前文設置計算至穩(wěn)定燃燒階段(作為無激勵工況進行對照)，而后以此結果為初始狀態(tài)并開啟聲學激勵繼續(xù)仿真。根據(jù)2.1節(jié)無激勵設計工況下矩形燃燒室構型的計算結果可知，燃燒室內的混合比O/F分布極不均勻，主反應區(qū)和高溫燃氣大部分聚集在燃燒室前半段，而聲學激勵的存在將進一步促進甲烷與氧氣之間的反應，因此本文選取燃燒室前半段施加如圖2所示的在y方向上的速度擾動(見式(1))，其產生的影響相比于選取燃燒室后半段將更明顯。數(shù)值計算中通過加載自定義函數(shù)(User Defined Function, UDF)將動量源項Mom′(即網(wǎng)格單元內的動量變化率，見式(2))導入動量守恒方程來實現(xiàn)。

圖2 橫向速度擾動設置Fig.2 Setup of transverse velocity disturbance

(1)

(2)

其中，v′0為擾動幅值(60 m/s)，h為燃燒室高度(141 mm)，f為擾動頻率，ρ為氣體密度。本文重點研究相同擾動幅值下噴霧燃燒對不同擾動頻率(3 000 Hz、4 538 Hz和6 000 Hz)的動態(tài)響應，其中4 538 Hz是該燃燒室的一階橫向(1T)振蕩模態(tài)固有頻率(fc=c/(2h)，c為熱力計算得到的相應燃燒室聲速(1 280 m/s))。

2 計算結果及討論

2.1 無激勵工況

由于目前對LOX/GCH4針栓式發(fā)動機的燃燒流場缺乏了解，本節(jié)首先采用開窗熱試觀測結果對仿真計算進行驗證。無激勵驗證工況的噴注器環(huán)縫寬度h和噴孔直徑d分別為2.23 mm和0.56 mm，對應GCH4和LOX的流量分別為37.2 g/s和84.5 g/s，其余發(fā)動機結構參數(shù)及仿真設置與本文1.1節(jié)一致，開窗熱試采用文獻[14]提出的單相機方案對噴霧(背景光成像)和火焰(自發(fā)輻射成像)同步觀測。穩(wěn)定燃燒階段噴嘴附近的噴霧軌跡和火焰形態(tài)如圖3所示，對比表明：①仿真得到的液氧液滴分布區(qū)域與熱試結果相近，均在針栓頭附近完成了快速蒸發(fā)，而觀察窗范圍內均以噴霧外緣區(qū)域燃燒最為劇烈，驗證了計算方法和網(wǎng)格的合理性；②熱試結果顯示液氧射流的噴霧錐角更大且射流表面存在表面波結構，原因在于仿真采用RANS方法和DPM分別模擬湍流流動和液氧射流，前者無法精確模擬連續(xù)相流場中的小尺度渦結構，后者忽略了液氧射流的破碎霧化過程(假設液氧在噴孔出口已霧化成離散液滴)，故無法反映真實的氣液霧化過程；③自發(fā)輻射成像得到的火焰分布范圍比仿真結果更廣，原因除了噴霧錐角的差異外，還在于試驗拍攝得到的圖像是沿視圖方向的積分，所有液氧射流蒸發(fā)燃燒都會在圖像上產生疊加作用，另外觀察窗玻璃邊緣位置的燒蝕現(xiàn)象也對圖像產生了干擾(如圖3所示紅色虛線框)，而仿真結果給出的僅是燃燒室對稱面上的歸一化釋熱率，主要體現(xiàn)上下兩股液氧射流的燃燒反應，軸線附近沒有明顯的火焰；④熱試圖像中火焰附著在液氧噴孔位置，而仿真得到的燃燒火焰距液氧噴孔有一段距離，原因在于仿真過程中燃燒反應的快慢由氣相組元的混合速率決定，液氧離開噴孔后要經(jīng)歷蒸發(fā)與混合后才能參與燃燒。

圖3 關于噴霧軌跡和火焰形態(tài)的試驗與仿真結果對比Fig.3 Comparison of experimental and simulation results on spray trajectory and flame distribution

下面針對1.1節(jié)給定的無激勵設計工況進行分析，圖4所示為該工況的燃燒室(P3)和集氣腔(P01)壓力振蕩波形，可以看出：燃燒室壓力時均值為1.71 MPa，基本達到設計值(1.8 MPa)；燃燒室中段壁面處的無量綱壓力振蕩幅度P′(1.17%)強于其他位置，未出現(xiàn)燃燒不穩(wěn)定現(xiàn)象(P′>3%)；集氣腔壓力振蕩略滯后于燃燒室壓力，且無量綱振蕩幅度P′僅為0.15%，這說明燃燒室壓力振蕩前傳至集氣腔的過程中出現(xiàn)了極大的耗散。圖4中燃燒室(P3)壓力振蕩波形表現(xiàn)出37 Hz和73 Hz這兩個振蕩主頻，根據(jù)文獻[15]可知，這兩個振蕩主頻分別由火焰在y方向的擺動過程和液氧的非穩(wěn)態(tài)蒸發(fā)過程引起的。

圖4 無激勵工況下P1與P3處達到穩(wěn)定的壓力振蕩波形Fig.4 Stable pressure oscillation wave at P01 and P3 without acoustic excitation

為了區(qū)分燃燒室內的火焰模式，這里引入火焰因子FI的概念[19]，如式(3)所示。

(3)

從圖4選取相隔約半個噴霧火焰擺動周期的兩個時刻，相應燃燒流場(見圖5)中GCH4來流集中在針栓頭下游并維持一定長度的低溫區(qū)，該低溫區(qū)伴隨燃燒火焰在y方向往復擺動，伴隨上下壁面附近回流區(qū)尺寸的此消彼長。初始噴注粒徑較小(10 μm)的液氧射流在穿透GCH4來流的過程中發(fā)生偏轉并快速蒸發(fā)，未表現(xiàn)出明顯的擺動現(xiàn)象。在回流區(qū)作用下高溫GOX往低溫區(qū)擴散并與GCH4發(fā)生劇烈反應，導致燃燒室內的混合比O/F分布極不均勻。另外，火焰因子的分布也表明燃燒室內的釋熱反應以擴散火焰為主，且火焰分布區(qū)域主要集中于兩條混合比等值線(0.5和500)之間。

(a) 內部溫度云圖及綠色的混合比等值線分布(0.5和500)(a) Interior temperature distribution with mixing ratio iso-lines (0.5 and 500) colored by green

2.2 擾動頻率對壓力振蕩的影響

在t=370.35 ms時刻開啟聲學激勵，不同擾動頻率下燃燒室各監(jiān)測點的無量綱壓力振蕩幅值如圖6所示，當擾動頻率與燃燒室1T振蕩模態(tài)固有頻率接近時壓力振蕩最為劇烈。位于燃燒室頭部的監(jiān)測點(P1和P2)的壓力時間歷程見圖7，f=4 538 Hz時開啟速度擾動至穩(wěn)定振蕩階段(振蕩幅度并未完全保持不變)所需的時間間隔最長(約50 ms)，相應的壓力時均值漲幅最大；在t=475.35 ms時刻關閉速度擾動后燃燒室壓力振蕩逐漸消失且壓力時均值開始下降，約50 ms后壓力振蕩波形恢復至無激勵狀態(tài)，這說明通過添加速度擾動可產生燃燒室聲學振蕩響應。隨著速度擾動的開啟和關閉，集氣腔壓力出現(xiàn)了相似的變化過程(以f=4 538 Hz為例，見圖8)，但產生的振蕩幅度和時均值的變化幅度均遠小于燃燒室壓力，這說明燃燒室壓力振蕩前傳至集氣腔的過程中出現(xiàn)了極大的耗散，且燃燒室壓力振蕩并非GCH4噴注引起。

圖6 不同擾動頻率下各監(jiān)測點的無量綱壓力振蕩幅度Fig.6 Dimensionless pressure oscillation amplitude of each monitoring point for different disturbance frequency

(a) f=3 000 Hz

圖8 f=4 538 Hz激勵條件下集氣腔內的壓力時間歷程Fig.8 Time histories of pressure oscillation in the gas manifold under acoustic excitation with f=4 538 Hz

結合圖7中局部放大圖和圖9所示不同擾動頻率下穩(wěn)定振蕩階段監(jiān)測點P1處的壓力頻譜可知，擾動頻率變化并未改變燃燒室壓力的1T振蕩模態(tài)，振蕩頻率均以擾動頻率為主，其中f=4 538 Hz時振蕩主頻對應幅值最大并產生高階諧振。

2.3 激勵條件下的聲學響應特性分析

下面結合燃燒流場和聲能平衡方程來解釋不同擾動頻率下噴霧燃燒對聲學激勵動態(tài)響應的差異。從圖7中不同擾動頻率下的壓力振蕩波形放大圖中選取相隔半個壓力振蕩周期的兩個壓力極值點對應的時刻，得到如圖10所示的燃燒流場，其中左列包括縱向截面上的溫度云圖和速度場矢量(以藍色箭頭表示)，右列包括液滴分布和z=0平面上的壓力場云圖，可以看出：激勵條件下燃燒室上下壁面附近的相對壓力變化符合1T振蕩模態(tài)，同時速度擾動促進了液氧蒸發(fā)和燃燒室前半段GOX/GCH4之間的混合反應，低溫區(qū)長度縮短。相比于3 000 Hz和6 000 Hz，f=4 538 Hz時噴霧燃燒過程對橫向聲學激勵的響應最為強烈，燃燒室內的回流區(qū)結構變化最大，燃燒反應在前半段基本完成，因此溫度進一步趨于均勻。Hardi等[9]通過旋轉齒輪控制輔噴管出口的開啟與關閉，從而對LOX/GH2同軸噴注縮尺火箭發(fā)動機施加燃燒室橫向聲學激勵，發(fā)現(xiàn)當激勵頻率等于燃燒室1T振蕩模態(tài)固有頻率時，反應速率(反比于火焰長度)和火焰區(qū)域的氫氧基(OH*)自發(fā)輻射強度顯著強于非共振頻率，與本文仿真結果一致。

圖9 不同擾動頻率下P1處的壓力頻譜圖Fig.9 Spectrogram plots of the pressure trace at P1 for different disturbance frequency

圖10 不同擾動頻率下相隔半個1T模態(tài)壓力振蕩周期的燃燒流場Fig.10 Snapshots of combustion ow eld with an interval of about half a period of 1T-mode pressure oscillation for different disturbance frequency

另外，在速度擾動作用下燃燒火焰與壓縮波(通過速度矢量表示)在y方向發(fā)生同步擺動，其中以f=4 538 Hz更為明顯，根據(jù)Rayleigh準則，燃燒釋熱將為高壓區(qū)(波腹)提供能量并維持聲學振蕩模態(tài)。相比之下，針栓噴嘴下游的低溫區(qū)與液氧液滴分布則在y方向上僅出現(xiàn)同步的小幅擺動現(xiàn)象，原因在于液氧和GCH4來流分別在密度和軸向噴注速度上遠大于高溫燃氣(慣性不同)，聲學激勵對兩者的影響較小。

為了進一步定量分析壓力和燃燒釋熱之間的關系，在穩(wěn)定振蕩階段針對圖10中綠色虛線框(0≤x≤0.07 m，-0.07 m≤y≤0)內的流場區(qū)域進行積分求解面積平均壓力和釋熱率，壓力和釋熱率之間的關系如圖11所示，擾動頻率f取3 000 Hz、4 538 Hz和6 000 Hz時指定區(qū)域內的平均壓力和燃燒釋熱的同頻振蕩相位差分別為65.45°、-8.18°和-63.53°，相位差約等于(fc-f)·π/fc，因此在合適的擾動幅值下，相位差由擾動頻率與1T模態(tài)振蕩固有頻率之間的相對大小決定。

(a) f=3 000 Hz

(4)

式中，p為壓力，Q為燃燒釋熱，u為速度矢量，ρ為密度，“-”和“′”分別表示平均值和脈動值。

三種擾動頻率下壓力和燃燒釋熱的同頻振蕩相位差在±90°之間，這說明聲學激勵項S1作用下燃燒釋熱產生項S2隨時間的積分為正，聲學源項S大于聲能耗散項D，進而維持了燃燒室內的壓力振蕩幅度；在p′u′項的作用下激勵區(qū)往無激勵區(qū)傳遞聲能，因此燃燒室壓力振蕩幅度沿x方向顯著減小。當f=4 538 Hz時壓力和燃燒釋熱基本表現(xiàn)出同相同步的周期性振蕩，因此燃燒釋熱產生項S2及可維持的燃燒室壓力振蕩幅度均更大，同時該工況下燃燒反應最為充分，外加額外引入的動量源項，使得燃燒室壓力時均值超過了設計壓力(1.8 MPa)；而當關閉速度擾動后，由于燃燒釋熱產生項S2和流體流動產生項S3無法維持，在聲能耗散項D作用下燃燒室壓力迅速恢復至無激勵條件下的穩(wěn)定燃燒狀態(tài)，這說明該針栓式發(fā)動機對燃燒室聲學振蕩具有較強的耗散作用，需采取額外的措施(如本節(jié)引入的動量源項)才能維持噴霧燃燒與聲學振蕩之間的強耦合關系。

最后，以f=4 538 Hz為例分析聲學激勵對低頻振蕩的影響。如圖7(b)所示穩(wěn)定振蕩階段的壓力振蕩幅度的變化對應一個58 Hz的低頻振蕩(振幅過小，圖9中沒有體現(xiàn))，圖12給出了相隔約半個低頻振蕩周期的兩個時刻的燃燒流場，可以看出，低溫區(qū)在y方向的擺動幅度較大，而兩條混合比等值線(0.5和500)的分布區(qū)域相比于無激勵工況(見圖5)顯著減小，這進一步說明燃燒反應基本在燃燒室前段完成，并導致激勵條件下低頻振蕩對應的頻率(58 Hz)高于無激勵工況(37 Hz)，同時混合比等值線在速度擾動作用下出現(xiàn)小尺度波動。火焰因子的分布則顯示燃燒室內的釋熱反應表現(xiàn)出擴散燃燒和預混燃燒并重的特征。另外，火焰除了集中于針栓頭附近外，還在速度擾動的作用下廣泛分布于燃燒室的一側壁面附近。

圖12 f=4 538 Hz激勵條件下相隔約半個低頻振蕩周期的燃燒室內部溫度場(包含綠色混合比等值線(0.5和500))及對稱面上的燃燒因子分布Fig.12 Instantaneous interior temperature elds, superimposed on the contour are the mixing ratio iso-lines (0.5 and 500), and FI distributions in the combustor symmetry plane with an interval of about half a low-frequency oscillation period under acoustic excitation with f=4 538 Hz

3 結論

本文采用Euler-Lagrange方法對LOX/GCH4針栓式發(fā)動機在不同擾動頻率的橫向聲學激勵作用下的聲學振蕩開展數(shù)值模擬，并通過開窗熱試證明該計算方法和網(wǎng)格能較為準確地捕捉噴霧軌跡和火焰形態(tài)。主要結論如下：采用的聲學激勵措施能產生同頻率的燃燒室1T振蕩模態(tài)聲學響應，噴霧燃燒對聲學激勵的響應強弱受擾動頻率的影響較大，其與燃燒室1T振蕩模態(tài)固有頻率(4 538 Hz)的相對大小決定了壓力和燃燒釋熱相位差；相比無激勵工況，速度擾動使得液氧蒸發(fā)和GOX/GCH4之間的混合反應均有所增強，低溫區(qū)長度縮短，存在預混火焰轉變?yōu)閿U散火焰的趨勢，進而增大了低頻振蕩對應的頻率；速度擾動關閉后聲學源項無法維持導致燃燒室壓力迅速恢復至無激勵條件下的穩(wěn)定燃燒狀態(tài)，這說明該針栓式發(fā)動機對燃燒室聲學振蕩具有較強的耗散作用；由于聲能自動從激勵區(qū)往無激勵區(qū)傳遞，穩(wěn)定振蕩階段燃燒室內的壓力振蕩幅度沿x方向呈現(xiàn)顯著減小的趨勢。當擾動頻率等于4 538 Hz時達到穩(wěn)定振蕩階段所需時間間隔最長(約50 ms)，此時壓力和燃燒釋熱相位差僅-8.18°，兩者接近同相振蕩，可維持的無量綱壓力振蕩幅度高達0.36，且在聲學激勵這一動量源項作用下壓力時均值增長至無激勵工況的設計壓力(1.8 MPa)以上；同時回流區(qū)結構顯著改變，噴霧和火焰伴隨壓力橫向振蕩出現(xiàn)明顯的同步擺動現(xiàn)象。