馬列波，聶萬勝，馮 偉，豐松江

?

液氧煤油火箭發(fā)動(dòng)機(jī)不穩(wěn)定燃燒過程的數(shù)值分析

馬列波，聶萬勝，馮 偉，豐松江

（裝備學(xué)院航天裝備系，北京，101416）

采用歐拉-拉格朗日方法對液氧煤油發(fā)動(dòng)機(jī)燃燒室內(nèi)的兩相燃燒過程進(jìn)行數(shù)值模擬，在驗(yàn)證模型可靠性基礎(chǔ)上，分析無隔板工況下自激1階切向高頻不穩(wěn)定性燃燒出現(xiàn)的原因。結(jié)果表明：噴嘴間霧化錐發(fā)生相互干涉使得推進(jìn)劑空間分布不均，導(dǎo)致脈動(dòng)釋熱，同時(shí)燃燒室內(nèi)無隔板時(shí)橫向壓力波阻尼特性降低，使得燃燒室內(nèi)出現(xiàn)1階切向燃燒不穩(wěn)定性；在不穩(wěn)定燃燒過程中，壓力振蕩波形和頻率與釋熱波動(dòng)的波形和頻率產(chǎn)生耦合，耦合程度越高，所含釋熱波峰峰值數(shù)量越少，其振蕩幅值也將越大。

液氧煤油火箭發(fā)動(dòng)機(jī)；高頻燃燒不穩(wěn)定性；壓力振蕩；釋熱波動(dòng)

0 引 言

在美國F-1液氧煤油火箭發(fā)動(dòng)機(jī)研制過程中出現(xiàn)了噴注器面燒毀現(xiàn)象，其燒蝕痕跡有著明顯的徑向流，表明存在高頻不穩(wěn)定燃燒。為了解決該問題，共進(jìn)行了2 000多次的全尺寸熱試車[1]；文獻(xiàn)[2]、文獻(xiàn)[3]對20世紀(jì)90年代之前的研究成果進(jìn)行了總結(jié)。高頻不穩(wěn)定燃燒是燃燒室內(nèi)聲學(xué)過程與推進(jìn)劑的噴射、霧化、蒸發(fā)、混合和化學(xué)動(dòng)力學(xué)等燃燒子過程中的一個(gè)或多個(gè)耦合的結(jié)果[2]。由瑞利準(zhǔn)則[4]可知，燃燒不穩(wěn)定性的產(chǎn)生與燃燒釋熱波動(dòng)和燃燒室內(nèi)壓力振蕩耦合過程相關(guān)，而在目前的研究中，對釋熱與壓力振蕩變化過程的分析較少；Huang等人[5～7]分析了釋熱波動(dòng)與壓力之間的振蕩波形及相位關(guān)系對單噴嘴燃燒室內(nèi)的縱向燃燒不穩(wěn)定性的影響，但對全尺寸發(fā)動(dòng)機(jī)中的切向高頻燃燒不穩(wěn)定性與二者之間的關(guān)系缺乏足夠的認(rèn)識(shí)。

本文針對大推力液氧煤油火箭發(fā)動(dòng)機(jī)內(nèi)的切向高頻燃燒不穩(wěn)定性問題，在初邊值條件不施加任何擾動(dòng)的情況下，通過對推力室燃燒流場的三維全尺寸仿真，對燃燒室內(nèi)壓力變化和分布特點(diǎn)、燃燒釋熱波動(dòng)和壓力自激振蕩變化關(guān)系進(jìn)行分析。

1 數(shù)學(xué)物理模型及計(jì)算方法

1.1 數(shù)學(xué)物理模型

氣相控制方程采用歐拉坐標(biāo)下的N-S方程[8]，其通用形式為

假設(shè)混合氣體滿足理想氣體狀態(tài)方程，多組分狀態(tài)方程計(jì)算如下：

液相采用離散顆粒模型進(jìn)行描述，其在Lagrangian參考系下通過顆粒運(yùn)動(dòng)方程計(jì)算其軌跡。根據(jù)顆粒受力分析，得到顆粒運(yùn)動(dòng)方程為

采用的LOX/Kerosene單步總包反應(yīng)如下：

1.2 數(shù)值方法

采用有限體積方法對控制方程組進(jìn)行離散，兩相之間的質(zhì)量、動(dòng)量和能量交換在氣相控制方程中通過相應(yīng)的源項(xiàng)實(shí)現(xiàn)。采用求解非穩(wěn)態(tài)可壓縮流的PISO算法對湍流兩相燃燒流場進(jìn)行耦合仿真計(jì)算，對流項(xiàng)采用迎風(fēng)格式，擴(kuò)散項(xiàng)采用二階中心差分格式，在近壁區(qū)雷諾數(shù)較低處采用壁面函數(shù)法[14]。

2 網(wǎng)格及邊界條件

2.1 網(wǎng)格生成

針對大推力的液氧煤油火箭發(fā)動(dòng)機(jī)，對其燃燒室進(jìn)行了三維建模和網(wǎng)格繪制（網(wǎng)格1），其隔板噴嘴為一周六徑，長度為40 mm。將隔板噴嘴平移到與其他主噴嘴平齊得到不帶隔板噴嘴的燃燒室模型（網(wǎng)格2）。網(wǎng)格模型見圖1。

圖1 網(wǎng)格模型

2.2 邊界條件設(shè)置

采用質(zhì)量入口邊界條件和壓力出口邊界條件，燃燒室壓力試驗(yàn)值為18 MPa，噴嘴壓降為1.2 MPa。煤油的霧化過程，其旋流的霧化錐通過噴霧錐角設(shè)為40°，旋流分?jǐn)?shù)設(shè)為0.5來表示，其噴射速度設(shè)為18.1 m/s，粒子平均直徑設(shè)為50 μm，粒子分布假設(shè)為均勻分布。燃燒室壁面采用無滑移絕熱條件，采用網(wǎng)格1和網(wǎng)格2進(jìn)行的燃燒室燃燒流場仿真所對應(yīng)的算例分別為算例1和算例2。

3 結(jié)果與討論

3.1 計(jì)算模型驗(yàn)證

a）壓力

b）振蕩幅值

圖2 壓力及其振蕩幅值隨時(shí)間變化

由圖2可以看出，壓力穩(wěn)定后，其平均室壓為16.7 MPa，與試驗(yàn)值為18 MPa誤差在7%左右，且其最大壓力振動(dòng)幅值小于5%，屬于穩(wěn)定燃燒。圖3給出的煤油液滴運(yùn)動(dòng)軌跡明顯存在周向的旋流運(yùn)動(dòng)（其視角為煤油噴射的反方向），這與實(shí)際霧化情況一致。仿真所得的霧化錐結(jié)果和平均室壓與試驗(yàn)結(jié)果符合較好，驗(yàn)證了本文所采用模型的正確性。

圖3 煤油液滴運(yùn)動(dòng)軌跡

3.2 壓力振蕩過程分析

為了觀察能量釋放過程，對反應(yīng)劇烈區(qū)域相關(guān)參數(shù)進(jìn)行監(jiān)控，監(jiān)控點(diǎn)位于噴注器面板下游50 mm，且靠近燃燒室壁面處。壓力及振蕩幅值隨時(shí)間變化曲線如圖4所示。圖4a為算例2監(jiān)測點(diǎn)處壓力隨時(shí)間的變化曲線，圖4b為監(jiān)測點(diǎn)處11 ms后，振蕩幅值隨時(shí)間變化的曲線。

a）壓力

b）振蕩幅值

圖4 壓力及振蕩幅值隨時(shí)間變化曲線

從圖4中可以看出，去掉隔板噴嘴后，在發(fā)動(dòng)機(jī)的燃燒過程中，壓力出現(xiàn)了大幅振蕩，其最大振蕩幅值超過了燃燒室平均室壓的10%，即該發(fā)動(dòng)機(jī)發(fā)生了不穩(wěn)定燃燒，且在發(fā)生大幅振蕩時(shí)，壓力隨時(shí)間的變化形態(tài)類似于陡峭的激波形式，表明此時(shí)振蕩較為劇烈。

對仿真結(jié)果得到的壓力振蕩數(shù)據(jù)進(jìn)行頻譜分析，得到的結(jié)果如圖5所示。

圖5 壓力數(shù)據(jù)的頻譜分析

由圖5可以看出，高頻區(qū)壓力振蕩頻率為1 235 Hz、1 512 Hz、2 027 Hz和2 832 Hz，表明在燃燒室中存在高頻不穩(wěn)定燃燒。另外，通過理論計(jì)算，該發(fā)動(dòng)機(jī)燃燒室的一階切向頻率為1 588 Hz，一階切向與一階縱向組合振型頻率為2 044 Hz。因此可以判斷燃燒室內(nèi)同時(shí)存在一階切向高頻不穩(wěn)定燃燒和一階切向與一階縱向組合振型不穩(wěn)定燃燒。仿真所得振蕩主頻值小于理論計(jì)算值，其原因?yàn)椋核捎玫幕瘜W(xué)反應(yīng)機(jī)理為單步總包反應(yīng)，忽略了燃燒過程中其他小分子產(chǎn)物，使燃?xì)猱a(chǎn)物平均分子量偏大，導(dǎo)致燃?xì)饴曀倨?，而燃燒室圓柱段切向聲學(xué)特性與燃燒室半徑及其燃?xì)饴曀傧嚓P(guān)，由于仿真采用的是全尺寸燃燒室模型，因此燃?xì)饴曀俚钠蛯?dǎo)致仿真結(jié)果低于理論值。圖6為一個(gè)振蕩周期燃燒室內(nèi)橫向壓力的分布（橫截面距離噴注面板5 mm），壓力分布形態(tài)呈現(xiàn)一階切向分布，高壓區(qū)域集中在靠近壁面的很小范圍內(nèi)，并且其隨時(shí)間沿著周向順時(shí)針移動(dòng)，由此也可以判斷燃燒室內(nèi)存在一階切向高頻不穩(wěn)定燃燒。

a）=24.4 ms

b）=24.8 ms

圖6 振蕩周期內(nèi)燃燒室橫向壓力的分布

c）=25.2ms

續(xù)圖6

反應(yīng)率是指單位時(shí)間參加化學(xué)反應(yīng)的推進(jìn)劑質(zhì)量，而該反應(yīng)所釋放的熱量則是其釋熱率，因此可以使用反應(yīng)率來表征釋熱率。本文將用反應(yīng)率替代釋熱率進(jìn)行分析，研究壓力振蕩與能量釋放之間的耦合過程。圖7給出了監(jiān)測點(diǎn)處在20~30 ms之間壓力與反應(yīng)率隨時(shí)間的變化曲線及10 ms之后壓力和反應(yīng)率振蕩數(shù)據(jù)的頻譜分析結(jié)果。從圖7a中的對比發(fā)現(xiàn)，壓力與反應(yīng)率均隨時(shí)間無規(guī)律振蕩，壓力的陡峭脈動(dòng)間斷性出現(xiàn)，在劇烈脈動(dòng)之后，其振蕩幅值顯著減小，部分壓力陡峭脈動(dòng)時(shí)反應(yīng)率也出現(xiàn)陡峭脈動(dòng)，表明兩者波形存在部分耦合。由圖7b中的頻譜特性對比發(fā)現(xiàn)，壓力振蕩頻率和反應(yīng)振蕩頻率在主頻1 235 Hz和 1 512 Hz處實(shí)現(xiàn)了耦合，并且壓力振蕩頻率2 027 Hz和2 832 Hz與反應(yīng)率振蕩頻率2 000 Hz和2 853 Hz非常接近。

a）壓力與反應(yīng)率隨時(shí)間的變化曲線

b）頻譜分析結(jié)果

圖7 壓力與反應(yīng)率隨時(shí)間變化曲線及其頻譜分析

噴嘴霧化場分布仿真結(jié)果如圖8所示。從圖8中發(fā)現(xiàn)（其視角為煤油噴射方向），部分噴嘴與噴嘴之間的噴霧錐存在相互干涉，使煤油液滴在燃燒室內(nèi)橫向分布不均勻，由此造成釋熱和壓力的初始脈動(dòng)，并且由于沒有隔板的存在，不能對壓力波的橫向傳播進(jìn)行抑制且降低了激發(fā)不穩(wěn)定燃燒所需的能量[15]，使壓力的自激振蕩得以發(fā)展，最終產(chǎn)生不穩(wěn)定燃燒。

圖8 噴嘴霧化場分布

在不穩(wěn)定燃燒階段，壓力與反應(yīng)率因推進(jìn)劑不均勻分布而造成的振蕩過程中，當(dāng)壓力振蕩與反應(yīng)率振蕩的波形和頻率耦合時(shí)，此時(shí)燃燒釋熱過程與壓力振蕩相位相同，壓力振蕩獲得能量，出現(xiàn)陡峭脈動(dòng)，其脈動(dòng)幅值體現(xiàn)了振蕩的劇烈程度。觀察發(fā)現(xiàn)，在壓力振蕩與反應(yīng)率振蕩波形耦合時(shí)，壓力脈動(dòng)幅值較大，其對應(yīng)的反應(yīng)率振蕩的波峰僅有一個(gè)峰值，壓力脈動(dòng)幅值較小，其對應(yīng)的反應(yīng)率振蕩的波峰含有多個(gè)峰值，表明當(dāng)反應(yīng)率振蕩的波峰僅有一個(gè)峰值時(shí)，壓力振蕩所獲得的能量全部用于增強(qiáng)其單一振蕩頻率波動(dòng)，能量利用率高，使得壓力振蕩幅值較大，振蕩更加劇烈；而當(dāng)反應(yīng)率振蕩的波峰含有多個(gè)幅值時(shí)，釋熱提供給壓力振蕩的能量分散在多個(gè)頻率的壓力波中，能量利用率低，使得壓力振蕩幅值較小。由于壓力振蕩波形與頻率只和部分反應(yīng)率振蕩波形與頻率實(shí)現(xiàn)耦合，壓力振蕩僅從部分釋熱中獲得能量，導(dǎo)致壓力振蕩除了維持自身的振蕩外，沒有多余的能量用于增強(qiáng)振蕩強(qiáng)度，其最大幅值僅超過平均室壓的12%，沒有達(dá)到其極限振蕩幅值。

4 結(jié) 論

對液氧煤油液體火箭發(fā)動(dòng)機(jī)燃燒室內(nèi)三維非穩(wěn)態(tài)兩相湍流燃燒過程進(jìn)行了數(shù)值仿真，得到無隔板工況下自激的一階切向高頻不穩(wěn)定燃燒仿真結(jié)果，發(fā)現(xiàn)在燃燒過程中，噴嘴間霧化錐發(fā)生相互干涉使得推進(jìn)劑空間分布不均，導(dǎo)致釋熱和壓力脈動(dòng)，同時(shí)由于沒有隔板的存在，不能對壓力波的橫向傳播進(jìn)行抑制，降低了激發(fā)不穩(wěn)定燃燒所需的能量，使得燃燒室內(nèi)壓力的自激振蕩得以發(fā)展，產(chǎn)生不穩(wěn)定燃燒；在不穩(wěn)定燃燒過程中，壓力振蕩劇烈程度與其和釋熱波動(dòng)波形與頻率的耦合程度及其振蕩波峰含有的釋熱波峰峰值數(shù)量有關(guān)，耦合程度越高，所含釋熱波峰峰值數(shù)量越少，其振蕩幅值將越大，振蕩將更加劇烈。

[1] 孫敬良, 張寶炯. 液氧-煤油火箭發(fā)動(dòng)機(jī)的燃燒不穩(wěn)定性問題[J]. 推進(jìn)技術(shù), 1995, 16(3): 34-41.

[2] Yang V, William A E. Liquid rocket engine combustion instability[M]. New York: The American Institute of Aeronautics' and Astronautics,1995.

[3] Harrje D T. Reardon F H. 液體推進(jìn)劑火箭發(fā)動(dòng)機(jī)不穩(wěn)定燃燒[M]. 朱寧昌, 張寶炯, 譯. 北京: 國防工業(yè)出版社, 1980.

[4] Durox D, Schuller T, Noiray N, et al. Rayleigh criterion and acousticenergy balance in unconfined self-sustained oscillating flames[J]. Combustion and Flame, 2009, 156: 106-119.

[5] Huang C, Anderson W E, Harvazinski M E, et al. Analysis ofself-excited combustion instabilities using decompositiontechniques[R]. AIAA 2013-1007, 2013.

[6] Harvazinski M E, Xia G, Anderson W E, et al. Analysis ofself-excited combustion instability using a combination of twoandthree-dimensional simulations[R]. AIAA 2012-0782, 2012.

[7] Harvazinski M E, Anderson W E, Merkle C L. Combustioninstability diagnostics using the rayleigh index[R]. AIAA2011-5548, 201

[8] 王福軍. 計(jì)算流體力學(xué)分析[M]. 北京: 清華大學(xué)出版社, 2004.

[9] Kappei F, Lee J Y, Johnson C E, et al. Investigation of oscillatory combustion progresss in actively controlled[R]. AIAA-2000-3348, 2000.

[10] Feng S J, Nie W S. He B. Control effects of baffle on combustion instability in a LOX/GH2rocket engine[J]. Journal of Spacecraft and Rocket, 2010, 47(3): 419-426.

[11] Sazhin S S. Advanced models of fuel droplet heating and evaporation[J]. Progress in Energy and Combustion Science, 2006, 32(2): 162-214.

[12] Sirignano W A. Fluid dynamics and transport of droplets and sprays[J]. Joural of Fluids Engineering, 2000, 122(1): 189 -190.

[13] Libby P A, Williams F A. Turbulent reacting flows[M]. New York: Springer Verlag, 1980.

[14] 聶萬勝, 豐松江. 液體火箭發(fā)動(dòng)機(jī)燃燒動(dòng)力學(xué)模型與數(shù)值計(jì)算[M]. 北京: 國防工業(yè)出版社, 2011.

[15] 尕永婧, 張會(huì)強(qiáng), 王希麟. 隔板對燃燒室壓力高頻自激振蕩的抑制作用[J]. 清華大學(xué)學(xué)報(bào), 2012, 52(7): 1008-1012.

Numerical Analysis of the Instability Combustion in LOX/Kerosene Rocket Engine

Ma Lie-bo, Nie Wan-sheng, Feng Wei, Feng Song-jiang

(Department of Aerospace Equipment, Equipment Academy, Beijing, 101416)

Numerically simulated thetwo-phase reacting process in LOX/Kerosene Rocket Engine combustion chamber by using Eulerian-Lagrangian method, after verified the validity of the model, analyzed the reason why the self-triggered first-order tangential high frequency instability combustion show up when there is no baffle. The result indicates: the mutual interference of atomizing cone between injectors makes the unevenly distributed in space, that leads to pulsation heat release, meanwhile the damping characteristics is lower when there is no baffle in the chamber, leads to the first-order tangential high frequency instability combustion in the chamber; In the instability combustion, pressure oscillation waveform and frequency couple with heat release oscillation waveform and frequency, the higher the degree of coupling is, the less number of the peak heat release wave crest contain , the oscillation amplitude will be higher.

LOX/Kerosene rocket engine; High frequency instability; Pressure oscillation; Heat release oscillation

1004-7182(2017)06-0032-05

10.7654/j.issn.1004-7182.20170608

V43

A

2016-05-31；

2017-11-03

國家自然科學(xué)基金項(xiàng)目(51206185，91441123)

馬列波（1992-），男，碩士研究生，主要研究方向?yàn)橐后w火箭發(fā)動(dòng)機(jī)