汪鳳山， 姚兆普， 劉 陽, 張 榛， 蔡 坤， 王 平， 毛曉芳， 楊 柳

1. 北京控制工程研究所, 北京 100094 2. 北京市高效能及綠色宇航推進(jìn)工程技術(shù)研究中心, 北京 100094

0 引 言

甲基肼(MMH)是當(dāng)今廣泛使用的宇航控制推進(jìn)系統(tǒng)及武器推進(jìn)系統(tǒng)用的可貯存液體燃料，常與四氧化二氮(NTO)、發(fā)煙硝酸等組成雙組元液體推進(jìn)劑，應(yīng)用于各種運載火箭、戰(zhàn)略和戰(zhàn)術(shù)裝備中.其中，甲基肼和四氧化二氮推進(jìn)劑組合具有良好的自維持燃燒特性，很短的著火延遲，穩(wěn)定可靠的啟動點火特性[1]，被廣泛應(yīng)用于雙組元姿控液體火箭發(fā)動機.

由于雙組元姿控發(fā)動機主要功能是為衛(wèi)星或航天器在軌姿態(tài)控制提供所需的沖量，短脈沖工況是最主要的工作模式.發(fā)動機燃燒室內(nèi)部推進(jìn)劑化學(xué)反應(yīng)通常難以達(dá)到穩(wěn)定燃燒狀態(tài)，而且MMH/NTO毒性較強，難以通過實驗手段對其化學(xué)反應(yīng)燃燒過程進(jìn)行詳細(xì)分析和研究[2].因此，選擇合理的MMH/NTO化學(xué)反應(yīng)動力學(xué)模型，采用數(shù)值仿真方法對發(fā)動機燃燒室內(nèi)流場進(jìn)行仿真分析，最終結(jié)合發(fā)動機高模熱試車燃燒試驗對其進(jìn)行驗證，是一種對雙組元姿控發(fā)動機脈沖工作燃燒過程進(jìn)行研究的優(yōu)選方法，對揭示發(fā)動機燃燒室內(nèi)瞬態(tài)反應(yīng)流動機理有著重要的意義.目前，尚未見MMH/NTO發(fā)動機采用簡化反應(yīng)機理的脈沖工況燃燒過程仿真與高模試驗對比研究的文獻(xiàn)報道.

本文采用了MMH/NTO簡化化學(xué)反應(yīng)模型，采用數(shù)值方法對某型雙組元液體火箭發(fā)動機脈沖工作燃燒過程進(jìn)行了仿真分析，并結(jié)合高模熱試車試驗結(jié)果對仿真結(jié)果進(jìn)行了驗證.

1 數(shù)值仿真模型

1.1 氣相湍流反應(yīng)流N-S方程

對于MMH/NTO為基礎(chǔ)的多組分化學(xué)反應(yīng)混合物，其氣相反應(yīng)流N-S方程通用形式為

(1)

式中，φ是通用因變量，Γφ是輸運系數(shù)，Sφ是氣相本身的源項，SPφ是液相產(chǎn)生的源項.通過對方程進(jìn)行雷諾分解和Favre平均可以得到平均量控制方程，由于表征湍流脈動引起的動量、質(zhì)量及能量輸運的二階關(guān)聯(lián)項是未知的，湍流燃燒過程中的化學(xué)反應(yīng)速率也是未知的，因此需要通過兩相模型、湍流燃燒模型和化學(xué)反應(yīng)模型使方程組封閉.

1.2 液滴控制方程

假定推進(jìn)劑噴出噴嘴后立即霧化，噴霧液滴的初始位置就在噴嘴出口處.對于霧化過程忽略液滴的一次破碎和霧化過程[3-5]，直接按照經(jīng)驗關(guān)系式給出推進(jìn)劑液滴霧化尺寸的初始分布，即人為地設(shè)定推進(jìn)劑霧化效果.

選擇兩個考慮因素較多的公式來計算索泰爾平均直徑DSMD(單位μm)如下：

(2)

經(jīng)過計算，甲基肼液滴的索泰爾平均直徑為34 μm，四氧化二氮液滴的索泰爾平均直徑為25 μm.

噴霧液滴初始速度依賴于噴嘴的具體形式，采用近似方法，假定總速度與噴嘴壓降的關(guān)系為

計算得到液滴的初始速度見表1所示.

表1 推進(jìn)劑液滴噴出初始速度Tab.1 Initial droplet velocities of propellants

液滴運動方程如下所示：

(3)

式中，Cd為阻力系數(shù)，F(xiàn)x,Fy,Fz可以看作是虛質(zhì)量力作用，這個力用來加速環(huán)繞在液滴周圍的流體，其大小為：

在時間步長上進(jìn)行積分就可以得到顆粒軌道上每一點液滴的瞬時速度.當(dāng)Tvap≤Tp

(4)

其中kc為傳質(zhì)系數(shù),由經(jīng)驗公式計算

(5)

當(dāng)計算出Ni后，就可由下式計算液滴的質(zhì)量變化：

mp(t+Δt)=mp(t)-NiApMw,iΔt

(6)

當(dāng)液滴溫度達(dá)到沸點溫度Tbp，并且液滴質(zhì)量大于揮發(fā)組分質(zhì)量時，采用沸騰速率方程

(7)

在液滴溫度小于其蒸發(fā)溫度時，液滴經(jīng)歷惰性加熱階段，溫升方程如下：

(8)

在液滴溫度達(dá)到蒸發(fā)溫度時，液滴同時經(jīng)歷惰性加熱階段和蒸發(fā)過程，溫升方程如下:

(9)

式中，Tp為液滴溫度(K)，h為對流傳熱系數(shù)(W/m2·K)，T∞為連續(xù)相溫度(K)，hfg為汽化潛熱(J/Kg).

推進(jìn)劑液滴除直接蒸發(fā)外，還濺射到發(fā)動機燃燒內(nèi)壁面形成冷卻液膜，有效地降低燃燒室壁面溫度，減小壁面熱流密度，從而保護(hù)火箭發(fā)動機壁面在許用范圍之內(nèi)[6-9].

1.3 化學(xué)反應(yīng)燃燒模型

本文采用MMH/NTO化學(xué)反應(yīng)動力學(xué)模型Mech23[10-1]，包括23組分、20個基元反應(yīng)，是經(jīng)簡化后的多步化學(xué)反應(yīng)機理.

表2 Mech23反應(yīng)動力學(xué)模型Tab.2 Mech23 reaction kinetics model

1.4 幾何結(jié)構(gòu)與網(wǎng)格

考慮到推力室的二維軸對稱結(jié)構(gòu)特點，采用了如圖1所示的2D計算域，在近壁面對網(wǎng)格進(jìn)行了適當(dāng)加密，使Y+在1的量級，共有約35 000個網(wǎng)格.

圖1 推力室?guī)缀谓Y(jié)構(gòu)與網(wǎng)格Fig.1 Thrust chamber geometry and mesh

1.5 邊界條件與計算工況

由于雙組元液體火箭發(fā)動機采用擠壓式系統(tǒng)結(jié)構(gòu)向燃燒室內(nèi)噴注燃料和氧化劑.在發(fā)動機啟動瞬間，由于燃燒室尚未完全建壓，因此，推進(jìn)劑流量在開機瞬間變化劇烈.根據(jù)試驗測量結(jié)果，發(fā)動機氧燃質(zhì)量混合比保持1.65不變，瞬間推進(jìn)劑總流量變化曲線如圖2所示.

圖2 發(fā)動機開啟瞬間推進(jìn)劑流量Fig.2 Propellant transient flow rates for engine starting up

推力室入口采用質(zhì)量流動入口邊界條件，瞬態(tài)總流量由圖2中擬合公式確定,氧化劑和燃料質(zhì)量流量根據(jù)總流量和混合比確定.推力室壁面采用無滑移壁面邊界條件，液膜與壁面作用采用了Stanton和O’Rourke的壁面液膜模型.根據(jù)液滴的沸點、撞擊能量和壁面溫度區(qū)分4種不同作用方式：黏滯、鋪開附著、反彈、鋪開并濺射[12-13].

推進(jìn)劑霧化粒度根據(jù)經(jīng)驗公式進(jìn)行預(yù)估，燃料液滴的粒徑設(shè)置為34 μm，NTO液滴的粒徑設(shè)置為25 μm.仿真計算主要有4 ms、6 ms、8 ms、10 ms和16 ms 5個工況.

2 結(jié)果與討論

2.1 瞬態(tài)流動燃燒分析

圖3為推力室內(nèi)靜壓在不同時刻下的變化圖.從圖中可以看出，在0時刻，推力室內(nèi)部壓力恒為100 Pa，0～3.7 ms之間，燃燒室壓力逐漸升高，在3.7 ms達(dá)到最大值0.95 MPa.在3.7～16 ms之間，燃燒室壓力逐漸減小，并最終建立了穩(wěn)定的壓力，約為0.7 MPa.到16 ms計算結(jié)束時，燃燒室內(nèi)壓力都再無大的變化.在關(guān)機之后1.7 ms，燃燒室內(nèi)壓力平緩下降，直到17.7 ms左右計算結(jié)束時，下降到大約0.07 MPa.

圖3 脈沖工況下推力室壓分布圖Fig.3 Distribution of thrust chamber pressure in pulse mode

圖4 脈沖工況下推力室靜溫分布圖Fig.4 Distribution of thrust chamber static temperature in pulse mode

圖4為推力室內(nèi)燃?xì)忪o溫在不同時刻下的變化圖.從圖中可以看出，在0～2 ms之間，燃燒室內(nèi)熱量逐漸積聚，溫度不斷升高.2 ms時，燃燒室中軸線高溫區(qū)進(jìn)一步減小，近壁面處低溫區(qū)向后移動，這是由于燃?xì)鉁囟壬撸旱握舭l(fā)速率加快所致.2.5 ms時，中心高溫區(qū)消失，推進(jìn)劑液滴在壁面的反彈導(dǎo)致燃料和氧化劑蒸發(fā)混合，在剪切層內(nèi)燃燒，高溫區(qū)溫度繼續(xù)上升.3 ms時，當(dāng)燃燒室平均溫度達(dá)到一定程度，即2 000 K左右時，氣相化學(xué)反應(yīng)速度進(jìn)一步加快，溫度不斷升高，此時燃燒室近壁面處高溫區(qū)向喉部移動，在接近喉部處的近壁面處存在一層很薄的高溫區(qū).在3.5 ms時刻，近壁區(qū)高溫消失，燃燒室內(nèi)溫度分布已經(jīng)接近最終的狀態(tài)，此時推進(jìn)劑液滴的軌跡把高溫區(qū)分成了兩個部分，一部分在燃燒室中心位置，另一部分在推進(jìn)劑液滴軌跡近壁面一側(cè)，這個高溫區(qū)一直存在，在此后的3.5～ 16 ms之間，溫度場逐漸穩(wěn)定.在16 ms時刻切斷推進(jìn)劑供應(yīng)，燃燒室溫度開始下降，在16.6 ms時刻，高溫燃?xì)忾_始流出燃燒室，由于沒有了持續(xù)的反應(yīng)物來源，燃燒室頭部出現(xiàn)了低溫區(qū).由于燃燒室近壁面處的燃?xì)饬魉佥^慢，此處的高溫燃?xì)饬鞒鐾屏κ业乃俣容^慢，一直到17.7 ms左右，關(guān)機后的計算結(jié)束，燃燒室內(nèi)溫度大幅度下降.

2.2 推力室壁面液膜熱流密度變化

圖5給出了液膜在整個生存周期內(nèi)通過壁面液膜熱流總通量的變化曲線，從圖中可以看出，在0～5 ms時間段內(nèi)，液膜熱流通量很小，這是因為在這段時間內(nèi)液膜在燃燒室內(nèi)壁面內(nèi)分布面積很小，且液膜附近燃?xì)鉁囟炔桓?在5～8 ms時間段內(nèi)，液膜逐漸發(fā)展到整個燃燒室，液膜蒸發(fā)速率逐漸加快，近壁面氣相組分濃度大大增加，化學(xué)反應(yīng)速率也相應(yīng)加快，導(dǎo)致近壁面處溫度很快上升，液膜內(nèi)外兩側(cè)溫差變大，液膜熱流通量因此加大.在8～16 ms之間，液膜的質(zhì)量由于劇烈消耗，逐漸被蒸干，液膜分布范圍減小，導(dǎo)致了壁面液膜熱流總通量迅速減小.

圖5 液膜熱流總通量隨時間變化曲線Fig.5 Curve of liquid film thermal flux varying with time

從圖中可以看出，發(fā)動機在脈沖工況下，推力室壁面液膜在壁溫不高的情況下，處于高效的核態(tài)沸騰的狀態(tài)，隨著壁溫升高，液膜熱流通量增大；但當(dāng)壁溫升高到一定程度后，壁面液膜由核態(tài)沸騰轉(zhuǎn)變?yōu)槟B(tài)沸騰狀態(tài)，液膜與壁面之間熱流通量急劇下降，冷卻效果急劇惡化.因此，甲基肼/四氧化二氮發(fā)動機推力室壁面液膜冷卻效果受脈沖寬度影響較大，為保證發(fā)動機在軌應(yīng)用可靠性，應(yīng)根據(jù)實際脈沖工況對發(fā)動機液膜冷卻效果進(jìn)行分析和地面點火驗證.

2.3 不同脈沖寬度的室壓變化

從圖6可以看出，從2 ms開始，壓力增加幅度加大；在3.7 ms左右，室壓迅速攀升到最高點，接近0.95 MPa，然后回落調(diào)整在0.7 MPa上下脈動.在4,6,8,10和16 ms不同的脈沖寬度下，燃燒室內(nèi)的壓力上升過程相同，而且在各自的脈沖時間結(jié)束關(guān)機以后的壓力下降過程也基本相同，都要經(jīng)過1.8 ms左右的時間才能使得燃燒室壓力下降到原先穩(wěn)定水平的10%.但在4 ms脈沖寬度下，壓力還未回落到穩(wěn)定狀態(tài)即已關(guān)機.

圖6 燃燒室壓力隨脈沖寬度變化圖Fig.6 Curves of thrust chamber pressure varying with time

2.4 高模熱試車試驗驗證

通過將該型雙組元姿控發(fā)動機在42 km高模試車臺進(jìn)行脈沖工況熱試車試驗，對發(fā)動機脈沖工況下的沖量特性進(jìn)行考核驗證，試驗過程中發(fā)動試車照片如圖7示.

圖7 發(fā)動機脈沖工況點火照片F(xiàn)ig.7 Photo of engine firing in pulse mode

圖8為發(fā)動機在不同脈沖寬度下，仿真分析與高模熱試車測量值對比示意圖.

從圖中可以看出，仿真分析結(jié)果與試驗測量結(jié)果吻合得很好.在圖中可以看出，計算得到的沖量在3 ms左右出現(xiàn)拐點，斜率變大，此時推進(jìn)劑化學(xué)反應(yīng)基本穩(wěn)定.試驗測量得出的沖量曲線在2.5 ms左右出現(xiàn)斜率變化，表明在2.5 ms之后推力室內(nèi)化學(xué)反應(yīng)基本穩(wěn)定，能夠輸出較穩(wěn)定的推力.由于仿真分析中沒有考慮因推力室外壁面輻射冷卻導(dǎo)致燃?xì)鉁囟认陆档膿p失，仿真分析結(jié)果在4 ms后要略小于試驗測量值，在16 ms脈沖結(jié)束時計算沖量與試驗值偏差約15%，可在后續(xù)研究中加以優(yōu)化改進(jìn).

圖8 發(fā)動機脈沖工況仿真與試驗結(jié)果對比Fig.8 Comparisons of numerical simulation and test results

基于上述分析和驗證結(jié)果，針對本文研究的發(fā)動機，脈沖寬度在6～8 ms下，發(fā)動機沖量輸出穩(wěn)定，且壁面液膜冷卻效率較高，是獲取發(fā)動機穩(wěn)定的最小脈沖沖量單元的優(yōu)選工況.

3 結(jié) 論

本文采用仿真分析和試驗研究的方法研究了MMH/NTO發(fā)動機脈沖工況下流場和溫度場演進(jìn)過程，揭示了發(fā)動機脈沖啟動過程中液膜熱流密度和沖量隨時間的變化規(guī)律，并結(jié)合高模熱試車試驗結(jié)果對仿真結(jié)果進(jìn)行了對比驗證，二者吻合得很好.

研究結(jié)果表明，由于推進(jìn)劑著火延遲的影響和啟動瞬間推進(jìn)劑流量的不穩(wěn)定性，MMH/NTO發(fā)動機通常在閥門開啟3 ms之后才能達(dá)到相對穩(wěn)定的燃燒狀態(tài)，輸出較為穩(wěn)定的推力或沖量；另外由于推力室壁面液膜冷卻效果與壁面液膜冷卻效果與脈沖寬度和占空比關(guān)系較大，在應(yīng)用中需要根據(jù)需要對脈沖工況下發(fā)動機液膜冷卻效果進(jìn)行分析和地面點火驗證，并根據(jù)分析和驗證結(jié)果優(yōu)選發(fā)動機脈沖工況進(jìn)行優(yōu)選，以保證發(fā)動機在軌工作可靠性.