高壓氫氧火箭發(fā)動(dòng)機(jī)推力室燃燒穩(wěn)定性分析

劉 倩，李敬軒，孫紀(jì)國，梁炫燁，向小林，潘 亮，鄭孟偉

(1.北京航天動(dòng)力研究所,北京 100076； 2.北京航空航天大學(xué) 宇航學(xué)院，北京 100191)

0 引言

高壓氫氧火箭發(fā)動(dòng)機(jī)具有大推力、高比沖等系列優(yōu)點(diǎn)，是重型火箭芯級動(dòng)力的重要支撐。大推力高壓火箭發(fā)動(dòng)機(jī)單臺推力為百噸量級，燃燒室尺寸顯著增加，燃燒室發(fā)生高頻燃燒不穩(wěn)定的傾向大大增加，合理設(shè)計(jì)噴注器結(jié)構(gòu)對于抑制推力室燃燒不穩(wěn)定的發(fā)生具有重要意義。高頻燃燒不穩(wěn)定依據(jù)反饋機(jī)理通常分為兩類，固有機(jī)理燃燒不穩(wěn)定和噴注耦合燃燒不穩(wěn)定。在噴注耦合聲學(xué)振型中，推進(jìn)劑噴注壓力和流量振蕩起主要作用。而對于固有機(jī)理燃燒不穩(wěn)定，推進(jìn)劑噴注后的霧化、蒸發(fā)、釋熱等子過程振蕩起主要作用，通常需要采用隔板、聲腔等穩(wěn)定裝置加以抑制。

對于氫氧火箭發(fā)動(dòng)機(jī)，燃燒室穩(wěn)定裝置可以改進(jìn)不甚穩(wěn)定的設(shè)計(jì)方案的穩(wěn)定性，但對較穩(wěn)定的方案的穩(wěn)定性裕度基本沒有改進(jìn)。美國航天飛機(jī)SSME和日本LE-7補(bǔ)燃?xì)?氧發(fā)動(dòng)機(jī)設(shè)置隔板和聲腔僅為防患于未然和減小風(fēng)險(xiǎn)，實(shí)際上對無隔板和聲腔發(fā)動(dòng)機(jī)做脈沖鑒定，證明發(fā)動(dòng)機(jī)是動(dòng)態(tài)穩(wěn)定的。RD-0120、LE-7A高壓補(bǔ)燃?xì)?氧發(fā)動(dòng)機(jī)燃燒室沒有加隔板和聲腔，試驗(yàn)結(jié)果證明燃燒穩(wěn)定。對于氫氧補(bǔ)燃循環(huán)發(fā)動(dòng)機(jī)推力室，其噴前溫度遠(yuǎn)遠(yuǎn)超過不穩(wěn)定邊界溫度，一般不會出現(xiàn)燃燒不穩(wěn)定。

在氫氧高壓火箭發(fā)動(dòng)機(jī)噴注器設(shè)計(jì)中，通過優(yōu)選氧噴嘴的聲學(xué)頻率使其與燃燒室各振型聲學(xué)錯(cuò)頻，可有效避免噴注耦合不穩(wěn)定的發(fā)生。但是，在未加隔板和聲腔的現(xiàn)有噴注器結(jié)構(gòu)下是否會發(fā)生固有機(jī)理燃燒不穩(wěn)定尚未可知。因此，迫切需要開展氫氧高壓火箭發(fā)動(dòng)機(jī)噴注器固有機(jī)理燃燒穩(wěn)定性分析。

近年來，由于計(jì)算機(jī)性能、流體力學(xué)和燃燒學(xué)的高速發(fā)展，學(xué)者采用數(shù)值模擬開展燃燒穩(wěn)定性預(yù)測，復(fù)雜的反應(yīng)機(jī)理和燃燒室結(jié)構(gòu)對熱聲不穩(wěn)定的預(yù)測成本和工程應(yīng)用帶來極大挑戰(zhàn)。Urbano等采用可壓縮大渦模擬LES模擬分析了德國宇航中心的縮比BKD氫氧發(fā)動(dòng)機(jī)的燃燒不穩(wěn)定，每計(jì)算1 ms時(shí)間長度需要1×10核時(shí)的CPU時(shí)間。尕永婧等針對液氧/煤油火箭發(fā)動(dòng)機(jī)燃燒室開展了三維非穩(wěn)態(tài)兩相燃燒過程的數(shù)值仿真，研究表明霧化角為65°時(shí)易于在燃燒室中出現(xiàn)爆炸性的可燃?xì)鈭F(tuán)并導(dǎo)致劇烈的壓力振蕩?？紤]到數(shù)值仿真精度、時(shí)間成本和工程研制周期間的最優(yōu)匹配，具有快速預(yù)測燃燒穩(wěn)定性的數(shù)值方法在工程研制中極具應(yīng)用價(jià)值。不穩(wěn)定預(yù)測模型可以研究設(shè)計(jì)參數(shù)對燃燒不穩(wěn)定性的阻尼和激勵(lì)作用進(jìn)而預(yù)測發(fā)動(dòng)機(jī)的穩(wěn)定極限。文獻(xiàn)[18-19]采用解耦法成功預(yù)測了劍橋大學(xué)鈍體湍流燃燒實(shí)驗(yàn)臺的極限環(huán)振蕩燃燒不穩(wěn)定性、ORACLES長火焰燃燒室的燃燒不穩(wěn)定性，其均準(zhǔn)確預(yù)測了燃燒不穩(wěn)定現(xiàn)象，且頻率預(yù)測精度偏差均小于2%。該解耦方法能夠在保證預(yù)測精度的情況下，大大降低計(jì)算量。

本文充分借鑒該研究思路，以高壓氫氧火箭發(fā)動(dòng)機(jī)推力室作為研究對象，針對燃燒室熱試驗(yàn)中易出現(xiàn)的一縱一切等聲學(xué)固有頻率，開展非穩(wěn)態(tài)燃燒與聲學(xué)系統(tǒng)解耦的燃燒穩(wěn)定性預(yù)測仿真研究，采用解耦法研究高頻燃燒不穩(wěn)定激勵(lì)機(jī)理以及發(fā)展過程，分析不同噴注器結(jié)構(gòu)參數(shù)和工況參數(shù)下推力室燃燒穩(wěn)定性。所做工作將為高壓氫氧火箭發(fā)動(dòng)機(jī)噴注器設(shè)計(jì)及燃燒穩(wěn)定性裕度評估提供參考。

1 預(yù)測方法介紹與計(jì)算模型

1.1 燃燒穩(wěn)定性預(yù)測方法介紹

不穩(wěn)定燃燒的反饋回路如圖1所示，將非穩(wěn)態(tài)燃燒與聲學(xué)系統(tǒng)解耦，將聲學(xué)擾動(dòng)(表現(xiàn)為速度脈動(dòng)或當(dāng)量比脈動(dòng))下的燃燒響應(yīng)(一般為燃燒釋熱率脈動(dòng))表征為火焰?zhèn)鬟f函數(shù)。將得到的火焰?zhèn)鬟f函數(shù)代入聲學(xué)求解器，預(yù)測燃燒不穩(wěn)定性。燃燒穩(wěn)定性預(yù)測計(jì)算主要分為兩個(gè)部分：①采用非定常雷諾平均NS方程(unsteady Reynolds averaged Navier-Stokes equations，URANS)仿真燃燒過程來計(jì)算火焰?zhèn)鬟f函數(shù)；②采用COMSOL來計(jì)算燃燒室的聲學(xué)模態(tài)，結(jié)合計(jì)算出來的分布式火焰?zhèn)鬟f函數(shù)，來預(yù)測燃燒不穩(wěn)定性。

圖1 燃燒不穩(wěn)定反饋回路Fig.1 Combustion instability feedback loop

當(dāng)燃燒室內(nèi)的流動(dòng)速度相對較低時(shí)，可近似將燃燒與聲學(xué)計(jì)算解耦，此時(shí)對火焰的擾動(dòng)主要來自上游的流量擾動(dòng)(或速度擾動(dòng))。對于帶有短噴管的推力室，由于噴管內(nèi)流動(dòng)速度較大，為了滿足解耦條件，將噴管從馬赫數(shù)0.3的位置截?cái)啵⒊隹谠O(shè)為聲學(xué)阻抗邊界條件，通過Magnus展開法計(jì)算得到，噴嘴進(jìn)口可認(rèn)為是聲學(xué)閉口條件。由于噴嘴面積總和/燃燒室橫截面積遠(yuǎn)遠(yuǎn)小于1，且噴嘴內(nèi)與燃燒室的溫差較大，因此噴嘴上游的聲學(xué)系統(tǒng)可與燃燒室內(nèi)的聲學(xué)系統(tǒng)解耦。本模型的計(jì)算主體如圖2所示。

圖2 計(jì)算組成Fig.2 Computational composition

1.2 噴注單元燃燒過程數(shù)值模型

推力室采用液氧/富氫燃?xì)?氣氫作為推進(jìn)劑，采用同軸直流噴注單元同心圓式排列?？紤]到計(jì)算時(shí)間成本，并結(jié)合以往計(jì)算經(jīng)驗(yàn)，對單個(gè)噴注單元進(jìn)行URANS數(shù)值模擬。下游燃燒計(jì)算區(qū)域采用圓柱形，徑向尺寸采用等效直徑法，即其截面積等于所有噴嘴的火焰占據(jù)的總橫截面積除以噴嘴的數(shù)量，求得單個(gè)噴注單元燃燒區(qū)域的等效直徑?？紤]到軸向長度對燃燒流場的影響，并綜合考慮初步計(jì)算得到的火焰長度，燃燒區(qū)域長度選取了50倍的氧噴嘴通道直徑，見圖3。

圖3 單噴嘴計(jì)算域Fig.3 Single injector computational domain

通過穩(wěn)態(tài)計(jì)算獲得燃燒室燃燒流場，以此為初始狀態(tài)，加入質(zhì)量流量擾動(dòng)后對噴嘴熱釋放振蕩過程進(jìn)行瞬態(tài)計(jì)算?？紤]到實(shí)際燃燒室的進(jìn)口和出口聲學(xué)邊界條件非常復(fù)雜，壓力振蕩擾動(dòng)頻率采用最易出現(xiàn)的一縱一切模態(tài)預(yù)測頻率。

采用ICEM劃分非結(jié)構(gòu)網(wǎng)格，并進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證，計(jì)算了網(wǎng)格數(shù)量為1.04×10、2.03×10、4.30×10的3種網(wǎng)格尺度計(jì)算，對比燃燒室軸線上的密度分布曲線，其中軸向起始位置為氧噴嘴入口。發(fā)現(xiàn)2.03×10、4.30×10的兩種網(wǎng)格尺度計(jì)算結(jié)果接近，故后續(xù)選取2.03×10的網(wǎng)格進(jìn)行計(jì)算，見圖4。利用商業(yè)軟件Fluent，采用分離求解器、SIMPLEC算法對壓力與速度進(jìn)行耦合。采用-SST湍流模型、渦耗散模型(eddy-dissipation model，EDM)、化學(xué)反應(yīng)模型。超臨界氫、氧以及燃?xì)饩催B續(xù)相處理，其密度等物性參數(shù)采用SRK狀態(tài)方程，其余物性參數(shù)選取NIST的數(shù)據(jù)進(jìn)行溫度的分段多項(xiàng)式擬合。瞬態(tài)仿真時(shí)間步長取1×10s，滿足Courant-Fredreichs-Lewy(CFL)條件。

圖4 網(wǎng)格無關(guān)性驗(yàn)證Fig.4 Grid independence verification

1.3 火焰?zhèn)鬟f函數(shù)計(jì)算

由于液氧的密度要遠(yuǎn)遠(yuǎn)大于燃?xì)獾拿芏龋趼返淖杩挂策h(yuǎn)遠(yuǎn)大于燃?xì)饴返淖杩?，因此在發(fā)生不穩(wěn)定燃燒時(shí)，燃料路的流量振蕩遠(yuǎn)遠(yuǎn)大于氧路的流量振蕩，進(jìn)而引起火焰燃燒釋熱率發(fā)生振蕩，與燃燒室聲學(xué)系統(tǒng)耦合，產(chǎn)生燃燒不穩(wěn)定性。故在燃?xì)馔ǖ朗┘淤|(zhì)量流量擾動(dòng)，分析下游火焰流場各參數(shù)的變化，并計(jì)算火焰?zhèn)鬟f函數(shù)，開展燃燒不穩(wěn)定性預(yù)測。

在燃?xì)馔ǖ缹α髁渴┘诱覕_動(dòng)，入口擾動(dòng)形式為

(1)

對于液體火箭發(fā)動(dòng)機(jī)同軸直流噴嘴，其穩(wěn)態(tài)火焰長度與燃燒室內(nèi)聲波波長處于同一個(gè)量級，因此在預(yù)測燃燒穩(wěn)定性時(shí)火焰不能被認(rèn)為是軸向緊湊型的，需要考慮分布式火焰?zhèn)鬟f函數(shù)。通常當(dāng)分段后火焰長度不大于聲波波長的0.1，可認(rèn)為其遠(yuǎn)遠(yuǎn)小于聲波波長，符合緊湊火焰假設(shè)，可以開展燃燒穩(wěn)定性預(yù)測。

將火焰分段后各段的火焰?zhèn)鬟f函數(shù)定義為

(2)

(3)

(4)

從進(jìn)口質(zhì)量擾動(dòng)到參考點(diǎn)速度擾動(dòng)之間的傳遞函數(shù)為

(5)

結(jié)合式(4)和式(5)，即可獲得火焰?zhèn)鬟f函數(shù)式(2)。

1.4 燃燒室聲學(xué)模態(tài)計(jì)算

利用COMSOL軟件對整個(gè)燃燒室建立Helmholtz方程。通過耦合Helmholtz方程和火焰?zhèn)鬟f函數(shù)，結(jié)合進(jìn)、出口處的聲學(xué)邊界條件，求得整個(gè)系統(tǒng)的聲學(xué)模態(tài)，表達(dá)式為

(6)

將燃燒室?guī)缀误w與噴注單元六段式的分布式火焰的幾何模型進(jìn)行裝配，求解整個(gè)系統(tǒng)的Helmholtz方程，即可預(yù)估推力室內(nèi)各個(gè)模態(tài)的頻率和壓力分布。其中，計(jì)算區(qū)域介質(zhì)的平均溫度、平均密度和平均聲速采用URANS仿真得到的結(jié)果，比熱容比采用REFPROP軟件來計(jì)算。對幾何模型采用非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分，網(wǎng)格數(shù)為1.06×10，如圖5所示。

圖5 聲學(xué)模態(tài)計(jì)算模型及網(wǎng)格劃分Fig.5 Acoustic modal calculation model and meshing

燃燒室的固有頻率主要取決于燃燒室的結(jié)構(gòu)尺寸及其內(nèi)部介質(zhì)的熱力學(xué)參數(shù)，火焰擾動(dòng)對聲學(xué)模態(tài)頻率的影響較小，主要影響聲學(xué)模態(tài)的穩(wěn)定性。計(jì)算聲學(xué)模態(tài)的特征值為復(fù)數(shù)，實(shí)部為模態(tài)頻率，虛部為增長率。當(dāng)虛部大于0時(shí)，主要物理量如壓力脈動(dòng)幅值逐漸減小，燃燒室在該聲學(xué)模態(tài)上是穩(wěn)定的；當(dāng)虛部小于0時(shí)，壓力脈動(dòng)幅值逐漸增大，燃燒室在該聲學(xué)模態(tài)上是不穩(wěn)定的。

2 結(jié)果與討論

2.1 穩(wěn)態(tài)計(jì)算結(jié)果

通過穩(wěn)態(tài)燃燒計(jì)算得到火焰長度約為0.29 m，絕熱火焰溫度約為3 742 K，見圖6。穩(wěn)態(tài)計(jì)算流場參數(shù)分布曲線見圖7，其中軸向長度起始位置為燃燒室噴注面板。平均溫度沿軸線逐漸增加，到軸向長度220 mm后溫度基本保持不變；平均密度沿軸線先下降后趨于平緩；平均聲速呈先減小后增加趨勢，主要與組分和溫度等相關(guān)。初步預(yù)測的一縱一切模態(tài)的頻率為3 767 Hz(見圖8)，將該一縱一切模態(tài)計(jì)算頻率作為計(jì)算火焰?zhèn)鬟f函數(shù)時(shí)入口擾動(dòng)頻率的參考。

圖6 噴嘴燃燒溫度場分布Fig.6 Temperature distribution of injector combustion

圖7 單噴嘴燃燒室穩(wěn)態(tài)燃燒場計(jì)算Fig.7 Calculation of steady-state combustion field parameters of single-injector combustor

圖8 一縱一切模態(tài)壓力分布Fig.8 First-order longitudinal first-order tangential modal pressure distribution

2.2 非穩(wěn)態(tài)燃燒計(jì)算結(jié)果

根據(jù)燃燒室的平均聲速(1 562.1 m/s)和預(yù)測的一縱一切聲學(xué)頻率(3 767 Hz)，可得燃燒室聲波波長約為0.41 m。結(jié)合仿真計(jì)算火焰長度0.29 m，將火焰沿著軸向分為6段，利用此6段區(qū)域的熱釋放速率可近似代表燃燒室內(nèi)所有的熱釋放速率。

通過開展單噴注單元的非穩(wěn)態(tài)燃燒仿真計(jì)算，獲得不同參數(shù)在一個(gè)周期內(nèi)的分布云圖，見圖9。燃燒釋熱率主要分布在噴嘴出口處，成蝌蚪狀，并在約220 mm處消失。在一個(gè)周期內(nèi)，“蝌蚪頭部”的釋熱率幅值產(chǎn)生周期性變化，“蝌蚪尾巴”也在不斷擺動(dòng)。

圖9 一個(gè)周期內(nèi)x-y平面燃燒釋熱率的分布云圖Fig.9 Distribution cloud map of combustion heat release rate in x-y plane in one cycle

不同軸向位置6個(gè)區(qū)域的燃燒釋熱率脈動(dòng)見圖10。釋熱率振蕩曲線接近正弦曲線，但具有不同的振蕩幅值和相位，并且第一段的振蕩幅度最大，越往后振蕩幅度越小。

圖10 6段燃燒釋熱率脈動(dòng)幅值Fig.10 Six-stage combustion heat release rate pulsation

圖11為傳遞函數(shù)的幅值和相位，可以看出其幅值在第一段區(qū)域最大，然后逐漸降低。

圖11 6段質(zhì)量擾動(dòng)和釋熱率之間的傳遞函數(shù)Fig.11 Transfer function from mass perturbation to heat release rate for six regions

圖12為不同軸向位置從進(jìn)口質(zhì)量擾動(dòng)到參考點(diǎn)速度擾動(dòng)之間的傳遞函數(shù)的幅值和相位，其相位的變化較小。其中參考位置選為=-12 mm到=-4 mm的平均值。結(jié)合式(2)、式(4)、式(5)，獲得分布式火焰?zhèn)鬟f函數(shù)，如圖13所示。

圖12 不同軸向位置速度傳遞函數(shù)Fig.12 Velocity transfer functions at different axial position

圖13 各段火焰?zhèn)鬟f函數(shù)幅值相位圖Fig.13 Amplitude phase diagram of flame transfer function of each segment

2.3 燃燒穩(wěn)定性仿真結(jié)果

利用非穩(wěn)態(tài)燃燒與聲學(xué)系統(tǒng)解耦法計(jì)算加載了分布式火焰幾何模型，得到一縱一切模態(tài)壓力分布，如圖14所示。

圖14 帶分布式火焰?zhèn)鬟f函數(shù)的一縱一切模態(tài)壓力分布Fig.14 First-order longitudinal first-order tangential modal pressure distribution with flame transfer function

隨著推力室燃?xì)?氧噴注速度比的增加，燃燒室一縱一切模態(tài)頻率呈上升趨勢，模態(tài)增長率也總體呈上升趨勢，燃燒過程變得更加穩(wěn)定，見圖15。

圖15 噴注速度比對模態(tài)及其增長率的影響Fig.15 Effect of injection velocity ratios on mode and its growth rate

分析原因?yàn)殡S著燃?xì)?氧噴注速度比的增加，燃?xì)獾膰娮⑺俣仍黾?，噴嘴出口氣液相互作用增?qiáng)，液柱長度減短(見圖16)。

圖16 不同噴注速度比下的溫度分布云圖Fig.16 Temperature distribution under different injection velocity ratios

同時(shí)由于燃?xì)獾膰娮⑺俣仍黾?，燃?xì)庠谌紵抑械耐Ａ魰r(shí)間減短，燃燒室整體釋熱區(qū)后移，分布式火焰上最大釋熱率相應(yīng)減少(見圖17)，系統(tǒng)相對而言逐漸變得穩(wěn)定。

圖17 不同噴注速度比下的火焰?zhèn)鬟f函數(shù)Fig.17 Flame transfer function at different injection velocity ratios

隨混合比的增大，燃燒室一縱一切模態(tài)的頻率變化不大，但模態(tài)增長率逐漸降低，燃燒室的穩(wěn)定性變差，見圖18。

圖18 混合比對模態(tài)及其增長率的影響Fig.18 Effect of mixing ratios on mode and its growth rate

分析原因?yàn)樵诓煌旌媳认氯紵覂?nèi)的燃燒均為富燃狀態(tài)，隨著混合比的增加，燃燒室中的氧量相對增加，可消耗更多的燃料放熱，燃燒室內(nèi)溫度更高，即熱釋放變得更大，見圖19，在液體火箭發(fā)動(dòng)機(jī)燃燒室中，由于噴嘴附近為壓力波腹，因此此處釋熱強(qiáng)度增強(qiáng)將使得熱聲不穩(wěn)定系統(tǒng)變得更不穩(wěn)定。

圖19 不同混合比下的火焰?zhèn)鬟f函數(shù)Fig.19 Flame transfer function at different mixing ratios

燃燒室相對氧噴嘴壓降增大時(shí)，推力室一縱一切模態(tài)的頻率變化不大，增長率略有降低，但影響不大，見圖20。這是由于燃燒室火焰較長，為火焰敏感型不穩(wěn)定，壓降的改變主要為邊界條件的改變，對燃燒室穩(wěn)定性影響不大。

圖20 相對噴嘴壓降對模態(tài)及其增長率的影響Fig.20 Effect of relative nozzle pressure drop onmode and its growth rate

氧噴嘴縮進(jìn)深度比對推力室一縱一切模態(tài)的頻率基本無影響。當(dāng)縮進(jìn)深度比從0.42增加為0.67時(shí)，系統(tǒng)更趨于穩(wěn)定；但縮進(jìn)深度比繼續(xù)增加對穩(wěn)定性的影響效果相反，即存在一個(gè)最優(yōu)縮進(jìn)深度比，見圖21。

圖21 縮進(jìn)深度比對模態(tài)及其增長率的影響Fig.21 Effect of retraction depth ratio on mode and its growth rate

當(dāng)富氫燃?xì)鈬娗皽囟壬邥r(shí)，推力室一縱一切模態(tài)的頻率變化不大，模態(tài)增長率略有降低，但影響不大，見圖22。

圖22 噴前溫度對模態(tài)及其增長率的影響Fig.22 Effect of per-injection temperature on mode and its growth rate

3 結(jié)論

針對高壓氫氧火箭發(fā)動(dòng)機(jī)推力室不設(shè)置隔板噴嘴和聲腔的結(jié)構(gòu)方案，利用火焰?zhèn)鬟f函數(shù)+低階模型的解耦法，開展燃燒穩(wěn)定性裕度的仿真預(yù)測，針對燃燒室熱試驗(yàn)中易出現(xiàn)的一縱一切等聲學(xué)固有頻率，分析了不同噴注參數(shù)和結(jié)構(gòu)參數(shù)下燃燒室的燃燒穩(wěn)定性。主要結(jié)論如下：

1)在給定工況下預(yù)測得到的燃燒室均未出現(xiàn)燃燒不穩(wěn)定現(xiàn)象；

2)混合比對燃燒室一縱一切模態(tài)頻率影響較小，但隨混合比的增大，模態(tài)增長率逐漸降低，燃燒室的穩(wěn)定性變差；

3)相對氧噴嘴壓降對推力室一縱一切模態(tài)頻率及其模態(tài)增長率均影響較??；

4)氧噴嘴縮進(jìn)深度比對推力室一縱一切模態(tài)的頻率無影響，但存在一個(gè)最優(yōu)縮進(jìn)深度比使燃燒穩(wěn)定性裕度最高；

5)在推力室設(shè)計(jì)中通過增加燃?xì)?氧噴注速度比或降低燃燒室混合比，有利于提升燃燒穩(wěn)定性裕度。