雙組元推力器霧化燃燒過(guò)程數(shù)值模擬及液膜冷卻作用分析

曹 順，汪鳳山，陳 健

（北京控制工程研究所，北京 100190）

雙組元推力器霧化燃燒過(guò)程數(shù)值模擬及液膜冷卻作用分析

曹 順，汪鳳山，陳 健

（北京控制工程研究所，北京 100190）

在帶有液膜冷卻雙組元推力器的設(shè)計(jì)過(guò)程中，對(duì)燃燒室的數(shù)值仿真是一項(xiàng)非常重要的工作.充分考慮了常被忽略或簡(jiǎn)化處理的推進(jìn)劑霧化、液滴破碎、液膜形成等重要過(guò)程，選取EDC燃燒模型和rea lizable k-ε湍流模型，采用有限體積計(jì)算方法得到燃燒室內(nèi)部液滴分布、液膜分布、靜溫分布等重要數(shù)據(jù).最后，根據(jù)所得結(jié)果闡述了燃燒室內(nèi)的詳細(xì)工作過(guò)程，并著重分析了液膜的分布特點(diǎn)以及對(duì)壁面的冷卻作用，得到了兩個(gè)與液膜冷卻相關(guān)的結(jié)論.

雙組元；液膜冷卻；燃燒室；數(shù)值模擬

雙組元推力器是主要用于航天器的姿態(tài)控制、軌道控制、姿態(tài)機(jī)動(dòng)和位置保持等功能的執(zhí)行機(jī)構(gòu)，具有非常重要的作用.為了設(shè)計(jì)性能優(yōu)化的推力器，需要對(duì)燃燒室內(nèi)的霧化、液膜、燃燒、流動(dòng)等過(guò)程有清楚的認(rèn)識(shí).但燃燒室內(nèi)部工作過(guò)程極其復(fù)雜，工況非常惡劣，目前很難完全通過(guò)試驗(yàn)的方法得到溫度場(chǎng)、壓力場(chǎng)、速度分布等數(shù)據(jù).隨著計(jì)算流體力學(xué)、計(jì)算燃燒學(xué)等數(shù)值方法的逐漸成熟，仿真已經(jīng)成為雙組元推力器設(shè)計(jì)過(guò)程必要的一個(gè)環(huán)節(jié).

然而到目前為止，國(guó)內(nèi)對(duì)燃燒室的仿真絕大多數(shù)都忽略了液滴破碎、液膜形成等重要物理過(guò)程.總裝備部指揮學(xué)院的聶萬(wàn)勝雖然對(duì)雙組元推力器燃燒室內(nèi)做了大量的仿真工作［1-3］，但均沒(méi)有考慮推進(jìn)劑霧化過(guò)程，同時(shí)也忽略了液膜對(duì)燃燒室的影響；許坤梅博士對(duì)雙組元推力器的燃燒、流動(dòng)過(guò)程進(jìn)行仿真［4］，但沒(méi)涉及液滴二次破碎和液膜形成等重要過(guò)程；張鋒等［5］雖然分析了液膜對(duì)燃燒室壁面的冷卻作用，但是對(duì)推進(jìn)劑由液滴到液膜形成過(guò)程缺乏足夠的研究，也沒(méi)有得到液膜在壁面分布的詳細(xì)情況.本文總結(jié)了前人的工作經(jīng)驗(yàn)，充分考慮了推進(jìn)劑噴注、液滴破碎和液膜形成等因素，建立合理的燃燒模型和湍流模型，對(duì)國(guó)內(nèi)新近研制的25N雙組元推力器燃燒室的工作過(guò)程進(jìn)行了數(shù)值仿真，著重分析了液膜冷卻對(duì)推力器工作過(guò)程的重要影響，為推力器的設(shè)計(jì)提供了重要的參考.

1 數(shù)學(xué)模型與仿真計(jì)算

1.1 推力器物理過(guò)程描述

雙組元推力器內(nèi)部工作原理極其復(fù)雜，涉及到多個(gè)復(fù)雜過(guò)程的耦合.圖1簡(jiǎn)單介紹了雙組元推力器的基本結(jié)構(gòu)和工作原理.在忽略壁面?zhèn)鳠岬那疤嵯?，本文將推力室?nèi)部工作過(guò)程合理地簡(jiǎn)化為以下幾步：

1）推進(jìn)劑由噴注器噴出，在各種力的作用下連續(xù)液體破碎為細(xì)小液滴；

2）液滴與燃燒室內(nèi)壁相互作用，形成液膜；

3）液滴在運(yùn)動(dòng)過(guò)程中與氣相作用，蒸發(fā)為氣體；

4）NTO（N2O4）遇高溫快速分解為NO2；

5）MMH（甲基肼）與NO2反應(yīng)產(chǎn)生高溫燃?xì)猓?/p>

6）高溫燃?xì)庥扇紵医?jīng)擴(kuò)張噴管?chē)姵?

圖1 雙組元推力器內(nèi)部工作原理示意圖Fig.1 The working principle scheme of bipropellant thruster

1.2 連續(xù)相控制方程

雙組元推力器內(nèi)部的流動(dòng)可以通過(guò)N-S方程來(lái)描述.由于推力器采用雙漩渦的噴注方式，燃燒室內(nèi)燃?xì)鈳в泻軓?qiáng)的旋流.因此，本文采用帶旋流修正的realizable k-ε湍流模型.燃燒模型采用工程上應(yīng)用比較成熟的EDC模型.離散方式采用二階迎風(fēng)格式.連續(xù)相控制方程包括連續(xù)方程、動(dòng)量方程（忽略體積力）、能量方程、組元輸運(yùn)方程（組分 k），分別如式（1）～（4）所示：

式中，ρ為密度，t為時(shí)間，xi、xj為位移變量，ui、uj為速度變量，p為壓力，h為混合物靜焓值，hj、hk分別為組元j和k的靜焓，T為溫度，λeff為有效傳熱系數(shù)，jki是組分k的擴(kuò)散流量，τij為粘性應(yīng)力，Yk是組分k的質(zhì)量分?jǐn)?shù)，D為組分Yk的擴(kuò)散系數(shù).上述方程中的源相來(lái)源于離散相（液滴）和氣相燃燒.其中質(zhì)量源相Sm是由液滴蒸發(fā)產(chǎn)生的；動(dòng)量源相SM，i等于在一個(gè)網(wǎng)格內(nèi)所有液滴動(dòng)量的損失值；Sh包括由燃燒產(chǎn)生的能量和液滴蒸發(fā)的能量交換兩項(xiàng)；SYk是化學(xué)反應(yīng)的反應(yīng)產(chǎn)物.可假設(shè)氣相為理想氣體，其狀態(tài)方程為

式中，Ru為通用氣體常數(shù)，Yi為i組元的質(zhì)量分?jǐn)?shù)，Mi為i組元的分子量.

MMH與NTO之間的化學(xué)反應(yīng)采用一步總包反應(yīng)模型.其化學(xué)反應(yīng)方程式為

由于N2O4在常溫下就非常容易離解為NO2，在燃燒室內(nèi)的高溫條件下，其自身離解速度將遠(yuǎn)遠(yuǎn)高于N2O4與MMH的反應(yīng)速度.因此可以認(rèn)為N2O4在由液態(tài)轉(zhuǎn)化為氣態(tài)的一瞬間便全部離解為NO2.所以，燃燒室內(nèi)實(shí)際發(fā)生的化學(xué)反應(yīng)可近似為

1.3 離散相模型

離散相模型主要包括推進(jìn)劑霧化、液滴蒸發(fā)、液膜模型等.考慮到雙組元推力器的結(jié)構(gòu)和噴注方式，本文采用壓力漩渦噴注模型，液滴的破碎過(guò)程采用TAB液滴破碎模型.

液膜的形成過(guò)程非常復(fù)雜，一般可以這樣描述［6］：推進(jìn)劑液滴由噴注器噴出以一定的速度撞到燃燒室內(nèi)壁，與壁面相互作用，部分液滴被粘滯，部分鋪開(kāi)附著，部分反彈，還有一部分液滴破碎并濺射，其中附在壁面的液滴形成壁面液膜，液膜在各種力的作用下在邊緣處與壁面發(fā)生分離.在液膜形成過(guò)程中，液膜表面與氣相進(jìn)行相互作用，包括液膜蒸發(fā)、對(duì)流換熱、受到氣相的拖拽力，液膜底部與固體壁面通過(guò)傳熱進(jìn)行能量交換.其原理參見(jiàn)圖2.

圖2 液膜形成機(jī)理示意圖Fig.2 Scheme of liquid film formation mechanism

本文綜合考慮了雙組元推力器雙漩渦噴注方式、燃燒室內(nèi)工作特點(diǎn)以及工程應(yīng)用成熟度等因素，采用了Stanton和O’Rourke的壁面液膜模型［7-9］.此液膜模型充分考慮了液滴與壁面相互作用區(qū)域內(nèi)的主要信息，如液滴速度、物性、溫度、壁面溫度、壁面粗糙度等.模型描述的液滴與壁面的相互作用方式與雙組元推力器液膜形成機(jī)理相符.因而，采取Stanton和O’Rourke的壁面液膜模型是合理的.

基于液滴的撞擊能量和液滴的沸點(diǎn)區(qū)分四種不同作用方式.其中液滴撞擊能量定義為：

ρ是液滴密度，Vr是液滴速度，D是液滴直徑，σ是液滴表面張力，δbl是邊界層厚度，H0是液膜厚度.

當(dāng)液滴撞擊能量小于16J且液滴溫度低于其沸點(diǎn)，液滴會(huì)粘滯到壁面；當(dāng)液滴溫度高于其沸點(diǎn)時(shí)，并以一特定的角度和速度與壁面相互作用，液滴會(huì)鋪開(kāi)附著于壁面；當(dāng)液滴溫度和壁面溫度均高于液滴沸點(diǎn)，同時(shí)撞擊能量小于臨界撞擊能量Ecr（根據(jù)相關(guān)文獻(xiàn)Ecr取為57.7J［7］）時(shí)，液滴將發(fā)生反彈，其反彈能量恢復(fù)系數(shù)為：

其中：θi為液滴對(duì)壁面的撞擊角度；當(dāng)液滴溫度與壁面溫度均高于液滴沸點(diǎn)，且撞擊能量大于臨界撞擊能量Ecr時(shí)，液滴與壁面相互作用后，將濺射出多個(gè)新的液滴.在液膜的邊緣處，當(dāng)液膜邊緣應(yīng)力超過(guò)壁面對(duì)液膜的粘滯力時(shí)，液壁將發(fā)生分離.二者受力關(guān)系非常復(fù)雜，主要取決于液膜邊緣的液體粘度、溫度、速度以及表面情況等當(dāng)?shù)貤l件.關(guān)于液膜控制方程等其他詳細(xì)信息參見(jiàn)文獻(xiàn)［7-9］.

1.4 計(jì)算幾何模型

本文的計(jì)算模型選用國(guó)內(nèi)新近研制的25N雙組元推力器.推力器采用MMH/NTO推進(jìn)劑，設(shè)計(jì)比沖為295N·s/kg，噴管膨脹比為150，壁面粗糙度為5μm.

本文采用了二維軸對(duì)稱(chēng)的幾何模型，計(jì)算網(wǎng)格總數(shù)為281340.推力器初始條件見(jiàn)表1.燃燒室壁面設(shè)定為絕熱無(wú)滑移的邊界條件；燃?xì)饨?jīng)過(guò)推力室喉部后成為超音速流，超音速流下游不影響上游，在確定出口邊界條件時(shí)只需將內(nèi)點(diǎn)值直接賦于出口邊界上的流場(chǎng)參數(shù)即可.

表1 推力器穩(wěn)態(tài)工作初始條件Tab.1 Initial condition of thrust in stable working state

2 計(jì)算結(jié)果與分析

經(jīng)過(guò)數(shù)值計(jì)算，本文得到了雙組元推力器穩(wěn)態(tài)工況下燃燒室內(nèi)重要物理參數(shù)的分布數(shù)據(jù).如圖3所示，液態(tài)推進(jìn)劑由噴注器噴出后，在各種力的作用下破碎成大小不一的液滴（最大液滴直徑為408μm，最小直徑為1.27μm）.液滴表面在燃燒室高溫的作用下不斷的蒸發(fā)，由MMH和NTO摩爾分?jǐn)?shù)和化學(xué)反應(yīng)圖可以看出，未反應(yīng)的氣態(tài)推進(jìn)劑主要集中在燃燒室靠近噴注器的一側(cè)，二者在這里通過(guò)充分的摻混，發(fā)生化學(xué)反應(yīng)，放出大量的熱.其反應(yīng)速率最快的位置出現(xiàn)在噴注器出口附近，達(dá)到了337mol/（m3·s）.

圖3 燃?xì)馐异F化與化學(xué)反應(yīng)Fig.3 Particle atomization and reaction in combustion chamber

圖4給出了在絕熱壁面條件下燃燒室靜溫分布情況.在靠近噴注器的前半部分為化學(xué)反應(yīng)不完全區(qū)域，同時(shí)由于液滴蒸發(fā)過(guò)程也需要吸收大量的熱，因而其溫度相對(duì)后部較低；后半部分到喉部以前屬于化學(xué)反應(yīng)完全區(qū)域，絕大多數(shù)化學(xué)能轉(zhuǎn)化為燃?xì)獾臒崮?，其溫度相?duì)較高，最高溫度達(dá)到2880K.從圖中也可以看出越是靠近燃燒室內(nèi)壁面的區(qū)域溫度越低，這是由于液膜的冷卻效應(yīng)引起的.

圖4中標(biāo)記了兩點(diǎn)，A點(diǎn)和B點(diǎn).A點(diǎn)位于燃燒壁面，B點(diǎn)與A的橫向距離相同，縱向位置處燃燒室內(nèi)壁與軸線(xiàn)之間，此處溫度達(dá)到2880K.從圖5中看出A點(diǎn)的溫度約為2290K.A、B點(diǎn)的溫差接近600K，這說(shuō)明液膜對(duì)于燃燒室壁面的熱防護(hù)作用非常明顯.此外，從圖6中還可以看出，燃燒室壁面溫度是沿著軸線(xiàn)方向逐漸升高的，在喉部位置溫度最高，超過(guò)2400K.

圖4 燃燒室靜溫分布情況Fig.4 Static temperature distribution in combustion chamber

圖5 燃燒室內(nèi)壁溫度（絕熱壁面條件）Fig.5 Inwall temperature distribution（thermal insutated inwall）

圖6 放大了燃燒室壁面局部區(qū)域的溫度分布，進(jìn)一步說(shuō)明了液膜的冷卻效果.從圖中可以看到，越靠近壁面的區(qū)域溫度越低，由于液膜的冷卻作用，在緊貼壁面處會(huì)形成一層很薄的低溫層.

壁面溫度分布具有很強(qiáng)的波動(dòng)性，主要是由液膜分布的不均勻性和間斷性所致.而液膜這種不均勻和間斷性是因?yàn)槿紵覂?nèi)的燃燒和流動(dòng)進(jìn)程非常復(fù)雜，使得液滴的初始條件（速度、方向）各不相同，其與壁面的相互作用方式極其不規(guī)律造成的.

如果以噴注位置為軸線(xiàn)零點(diǎn)，液膜在壁面的總體分布區(qū)域從軸線(xiàn)方向0.35mm一直延伸至45mm左右.本文給出了4.6mm到8.6mm之間燃燒室內(nèi)壁液膜厚度的分布圖，如圖7所示.從圖中可以看出，A、B、C三點(diǎn)是明顯的液膜間斷處，不同位置的液膜厚度最大值為18.2μm.正是由于液膜分布的這種不均勻性，造成燃燒室壁面溫度分布的波動(dòng).

圖6 壁面區(qū)域靜溫分布放大圖Fig.6 Static temperature distribution in zone adjacent inwall（zoomed picture）

圖7 燃燒室內(nèi)壁面液膜分布（局部）Fig.7 Distribution of liquid film on the combustion chamber inwall（local）

3 結(jié) 論

本文對(duì)雙組元推力器燃燒室內(nèi)的工作過(guò)程進(jìn)行了數(shù)值仿真，考慮了在以往仿真工作中常被忽略的液滴破碎及液膜形成等重要進(jìn)程，得到了液滴霧化、化學(xué)反應(yīng)、靜溫、液膜分布等重要結(jié)果.根據(jù)數(shù)值仿真的結(jié)果，闡述了雙組元推力器燃燒室內(nèi)詳細(xì)的工作過(guò)程，著重分析了液膜的分布特點(diǎn)以及對(duì)燃燒室的冷卻效果，最后得到了兩個(gè)關(guān)于液膜冷卻的結(jié)論：

1）在雙組元推力器的穩(wěn)態(tài)工作過(guò)程中，液膜冷卻的效果明顯，對(duì)燃燒室壁面具有非常好的保護(hù)作用；

2）液膜在燃燒室內(nèi)壁的分布是間斷的，厚度不均勻（0～18.2μm），但這并不影響液膜的冷卻效果.

［1］ 聶萬(wàn)勝，莊逢辰.自燃推進(jìn)劑火箭發(fā)動(dòng)機(jī)穩(wěn)態(tài)燃燒過(guò)程的數(shù)值模擬［J］.推進(jìn)技術(shù)，1998，19（5）：7-9 Nie W S，Zhuang F C.Numerical simulation of stable combustion process in earth storable bipropellant rocket engine［J］.Journal of Propulsion Technology，1998，19（5）：7-9

［2］ 聶萬(wàn)勝，莊逢辰.噴霧特性對(duì)液體火箭發(fā)動(dòng)機(jī)燃燒穩(wěn)定性的影響［J］.推進(jìn)技術(shù)，2000，21（3）：56-59 Nie W S，Zhuang F C.Effect of spray characteristic on the liquid rocket combustion stability［J］.Journal of Propulsion Technology，2000，21（3）：56-59

［3］ 聶萬(wàn)勝，豐松，解慶紛，莊逢辰.液體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定性主動(dòng)控制數(shù)值仿真［J］.工程熱物理學(xué)報(bào)，2008，29（7）：1253-1255 Nie W S，F(xiàn)eng S J，Xie Q F，Zhuang F C.Numerical simulation of liquid rocket combustion instability active control［J］.Journal of Engineering Thermophysics，2008，29（7）：1253-1255

［4］ 許坤梅.雙組元小推力液體火箭發(fā)動(dòng)機(jī)推力室工作過(guò)程數(shù)值模擬［D］.北京：北京航空航天大學(xué)，2005 Xu K M.Numerical simulation of the working process in thruster chanber of bipropellant thruster［D］.Beijing： Beihang University，2005

［5］ 張鋒，仲偉聰.膜冷卻推力器室傳熱計(jì)算研究［J］.火箭推進(jìn)，2009，35（4）：35-37 Zhang F，Zhong W F.Computational investigation of heat transfer for film cooling thrust chamber［J］.Journal of Rocket Propulsion，2009，35（4）：35-37

［6］ Naber JD，Reitz R D.Modeling engine spray/wall impingement［J］.General Motors Research Laboratories，SAE Transactions，1989，97（6）：118-140

［7］ O’Rourke P J，Amsden A A.A particle numerical model for wall film dynamics in port-fuel injected engines［C］.International Fuels and Lubricants Meeting and Exposition，San Antonio，TX，Oct 14-17，1996

［8］ O’Rourke P J，Amsden A A.A spray/wall interaction submodel for the kiva-3 wall film model［R］.Detroit，SAE Technical Paper，2000-01-0271，2000

［9］ Stanton DW，Rutland C J.Modeling fuel film formation and wall interaction in diesel engines［R］.Detroit，SAE Technical Paper 960628，1996

Num erical Sim ulation of W orking Process in Cham ber of Bipropellant Thruster and Analysis of Film-Cooling Effect

CAO Shun，WANG Fengshan，CHEN Jian
（Beijing Institute of Control Engineering，Beijing 100190，China）

In the design of film-cooled bipropellant thruster，the numerical simulation of the combustion chamber is a necessary step.In this paper，process of atomizing of the liquid propellant，breaking up of drop lets，liquid film formation are taken into account，and the EDC combustion model and the rea lizable k-εturbulent flow model are selected.Distributions of droplets，liquid film and static temperature are obtained by using the numerical simulation.At last，based on the analysis of simulation results，the working process in combustion chamber are interpreted，and two conclusions on the liquid film cooling are derived.

bipropellant；liquid film cooling；combustion chamber；numerical simulation

V434

A

1674-1579（2012）02-0058-05

曹 順（1983—），男，碩士研究生，研究方向?yàn)楹教炱魍七M(jìn)技術(shù)；汪鳳山（1981—），男，工程師，研究方向?yàn)楹教炱魍七M(jìn)技術(shù)；陳 健（1969—），男，研究員，研究方向?yàn)楹教炱魍七M(jìn)技術(shù)。

2011-11-15

DO I：10.3969/j.issn.1674-1579.2012.02.012