兗立文　許坤梅　王慧潔

（北京航空航天大學(xué)宇航學(xué)院，北京　100191）

小推力液體火箭發(fā)動機(jī)沉降型液膜冷卻液膜長度研究

兗立文許坤梅王慧潔

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

沉降型液膜冷卻是指從噴嘴向壁面以一定角度噴注某種液滴，液滴到達(dá)壁面沉降形成液膜，實現(xiàn)對壁面熱防護(hù)的一種方法，常用于小推力液體火箭發(fā)動機(jī)熱防護(hù)系統(tǒng)中。采用Stechaman半經(jīng)驗方法對決定液膜冷卻效果的關(guān)鍵參數(shù)——液膜長度求解。采用k-w模型描述湍流流動、Eulerian-Lagrangian模型描述兩相流，采用C/C++語言編寫B(tài)ai-Gosman液滴撞壁模型程序，采用數(shù)值模擬推進(jìn)劑液滴撞壁后復(fù)雜的狀態(tài)變化過程，考慮壁面熱輻射，對400N雙組元MMH/NTO自然推進(jìn)劑發(fā)動機(jī)推力室內(nèi)的蒸發(fā)、流動、燃燒和傳熱過程進(jìn)行了數(shù)值模擬。在采用半經(jīng)驗公式方法求解液膜長度的過程中，分析了液膜流量對液膜長度的影響，研究了該方法的實際應(yīng)用情況。實驗結(jié)果表明，該方法可以較好地計算液膜長度，計算結(jié)果與實驗具有較高的一致性，對工程實踐具有重要的指導(dǎo)意義。

液體火箭發(fā)動機(jī)，液膜冷卻，液膜長度，液滴撞壁，數(shù)值模擬

引言

對于小推力液體火箭發(fā)動機(jī)來說，因為其推力和尺寸較小，所以，液態(tài)燃料噴入推力室后，不可避免地會出現(xiàn)液滴撞壁過程。液滴撞壁后，其運動過程非常復(fù)雜，在某些條件下會在壁面形成液膜。在對發(fā)動機(jī)推力室內(nèi)噴霧湍流燃燒進(jìn)行數(shù)值模擬的基礎(chǔ)上，采用Bai-Gosman［1］液滴撞壁模型模擬了液滴撞壁后的復(fù)雜狀態(tài)變化過程，并分析了液滴撞壁后的軌跡、推力室流場和發(fā)動機(jī)各項性能參數(shù)。

液體火箭發(fā)動機(jī)燃燒室溫度高、熱流密度大，可采用冷卻方法保護(hù)發(fā)動機(jī)，目前常用的冷卻方法有再生冷卻、輻射冷卻、熱沉法、隔熱法、液膜冷卻或發(fā)汗冷卻法［2］。其中，液膜冷卻是一種應(yīng)用較為廣泛且有效的冷卻方式。各國研究人員對液膜冷卻開展了多年的深入研究。例如，文獻(xiàn)［3～4］研究了推力室中心氣流與液膜界面的結(jié)構(gòu)，以及冷卻液膜的不穩(wěn)定性；文獻(xiàn)［5～6］從實驗和理論角度對液膜進(jìn)行了研究，但未對液膜長度進(jìn)行定量分析，且不針對小推力發(fā)動機(jī)。本文的研究對象為小推力發(fā)動機(jī)內(nèi)沉降型液膜冷卻，且需進(jìn)行定量分析。分析發(fā)現(xiàn)，Stechaman半經(jīng)驗方法［7］針對小推力發(fā)動機(jī)液膜冷卻進(jìn)行了研究，提出了適用于小推力發(fā)動機(jī)液膜長度計算的經(jīng)驗公式，因此，可采用該方法分析液膜長度，計算中需確定液膜流率，但不同于噴入型液膜冷卻，在沉降型液膜冷卻中，液膜流率是未知的。本文針對這一點，采用Bai-Gosman液滴撞壁模型，通過自編程求解液滴撞壁后形成液膜時的流率，采用Stechaman半經(jīng)驗方法對400N小推力發(fā)動機(jī)［8］的液膜冷卻過程進(jìn)行數(shù)值模擬，分析液膜流量因素對液膜長度的影響。

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

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

控制方程中的各種定律反映的都是單位時間單位體積內(nèi)物理量的守恒性質(zhì)，可以表示成以下通用形式［9］：

其中，Φ表示通用變量，可以代表u、v、w、T等求解變量；Г為廣義擴(kuò)散系數(shù)；S表示廣義源項。式（1）中各項依次為瞬態(tài)項、對流項、擴(kuò)散項和源項。

1.2液滴運動控制方程

液滴在空間的分布及其運動狀態(tài)在很大程度上取決于其運動過程中所受到的力。在忽略液滴的旋轉(zhuǎn)及流場中速度梯度產(chǎn)生的升力、Magnu力，以及重力等作用的情況下，液滴運動方程可以表示為：

通過對時間的積分，可得到液滴在各個時刻的速度和位置。CD是液滴阻力系數(shù)，rp是液滴的半徑，和分別是介質(zhì)氣體和液滴速度矢量。阻力系數(shù)CD由下式求得：

其中，Re為液滴的雷諾數(shù)。

1.3Bai-Gosman液滴撞壁模型

Bai和Gosman［2］提出的模型將液滴撞壁過程詳細(xì)劃分成粘附、反彈、鋪展、沸騰后破碎、反彈伴隨破碎、破碎、飛濺等7種形式。在以臨界韋伯?dāng)?shù)判定的濕壁情況下，不同狀態(tài)的轉(zhuǎn)變臨界值情況為：反彈→鋪展→飛濺。

式中，La-Laplace數(shù)據(jù)由壁面粗糙度決定，We為液滴的韋伯?dāng)?shù)，σ為液滴的彈性能，d為液滴的直徑，v為液滴的速度，ρ為液滴的密度，μ為液滴的粘度。本文僅討論液滴撞壁面后反彈和飛濺情況。

反彈后液滴的切向和法向速度如下：

其中，u和v分別是切向與法向分速度，下標(biāo)e和r分別代表入射和反彈，ε為恢復(fù)系數(shù)。ε=0.993-1.76θ+1.56θ2-0.49θ3，θ是以弧度計算的液滴入射角。

飛濺時，Bai-Gosman模型假定液滴撞壁后、飛濺后液滴團(tuán)的直徑和速度分別為d1、d2和U1、U2。如果兩個二次液滴團(tuán)所含的液滴數(shù)分別為N1和N2，其質(zhì)量守恒方程為：

借用氣象學(xué)中的相關(guān)數(shù)據(jù)擬合確定總的二次液滴數(shù)

二次液滴的速度方程為：

參考雨滴尺寸與速度的關(guān)系式推導(dǎo)出液滴的速度比為：

上式表明，若破碎后的液滴尺寸與速度呈負(fù)相關(guān)關(guān)系，尺寸越大，速度就越小。

液滴的動量方程為：

式中，Cf是摩擦系數(shù)，取值為0.6～0.8，θ1和θ2是相應(yīng)的角度。

1.4Stechaman半經(jīng)驗方法

在考慮推力室內(nèi)熱氣流與液膜的熱傳遞時，采用Stechaman［7］修正的Bartz［10］方法來計算熱氣流與液膜的傳熱系數(shù)［11］，傳熱系數(shù)為式（12）。

Stechaman半經(jīng)驗公式方法計算液膜長度為式（13）。

式中，*為喉部位置，μf為動力黏度，D為發(fā)動機(jī)燃燒室直徑，A為面積，W為冷卻劑流量，HW為壁面條件下的焓值，TW為壁面條件下的溫度，Hr為恢復(fù)焓值，Tr為恢復(fù)溫度，CPL為液膜冷卻劑的等壓比熱容，Prf為液膜的普朗特數(shù)，η為與液膜雷諾數(shù)相關(guān)的無量綱數(shù)，σ為邊界層處與密度和黏度有關(guān)的無量綱數(shù)，L為液膜長度，Tl為液滴初始噴注溫度，Ts為飽和溫度，P為發(fā)動機(jī)燃燒室周長，hg為熱氣流與液膜的傳遞系數(shù)，λ為冷卻劑的潛熱。

采用Stechaman半經(jīng)驗公式法計算液膜長度時，需確定形成液膜的流量，在噴入型液膜冷卻中，可根據(jù)入口條件獲得；而在沉降型液膜冷卻中，由于燃料從噴嘴向壁面以一定角度噴注某種液滴，而液滴撞壁后狀態(tài)復(fù)雜，不能直接確定其形成液膜的流量。因此，本文采用Bai-Gosman液滴撞壁模型，通過自編程求解液滴撞壁后形成液膜時的流率。

2　數(shù)值模擬

2.1計算參數(shù)

發(fā)動機(jī)的噴注器采用雙組元離心式，向其內(nèi)噴嘴注入燃料甲肼（MMH），外噴嘴噴入氧化劑硝基-三唑-酮（NTO），氧燃推進(jìn)劑流率的混合比為1.667，全部流量為130g/s，初始溫度為300K。液滴的其余參數(shù)由前面計算得出。MMH的噴注速度Vinj，MMH為16.2m/s，NTO的噴注速度Vinj，NTO為2.6m/s，與對稱軸的夾角為55°。液滴進(jìn)入燃燒室時，MMH的入口位置為：X=3.5mm，Y=5mm；NTO的入口位置為：X=3.5mm，Y=8mm。假定噴霧液滴的最小粒徑dmin為5μm，最大粒徑dmax為200μm或110μm；假設(shè)噴注器所產(chǎn)生的液滴尺寸分布服從Rosin-Rammler分布。按Rosin-Rammler分布將液滴分組，分組數(shù)目為50，則質(zhì)量分布函數(shù)公式為：

n為液滴分布指數(shù)，表示液滴直徑分布的均勻性，一般n=2～4。n越大，液滴直徑分布越均勻，本文取n=4。其液滴參數(shù)詳見表1。

表1　噴霧噴注參數(shù)

2.2計算網(wǎng)格

圖1是某400N軌控發(fā)動機(jī)的計算網(wǎng)格，共包括24600個網(wǎng)格單元。根據(jù)湍流對壁面附近網(wǎng)格的要求，對壁面附近的網(wǎng)格進(jìn)行合理加密。由于發(fā)動機(jī)燃燒室噴嘴入口附近流場的溫度、濃度等參數(shù)梯度較大，所以，需對噴注面附近的網(wǎng)格加密。

圖1　計算網(wǎng)格

2.3壁面邊界條件

氣相壁面為熱輻射無滑移邊界；液相壁面為Bai-Gosman液滴撞壁模型。

2.4數(shù)值解法

控制方程采用有限體積法進(jìn)行離散，離散方程的求解采用SIMPLE算法。圖2為計算方法流程：首先，對推力室內(nèi)流場進(jìn)行數(shù)值仿真，收斂時計算內(nèi)壁面溫度Tw， in［i-1］、獲取Stechaman經(jīng)驗公式法所需參數(shù)、計算液膜長度，更新設(shè)置壁面條件；其次，對推力室內(nèi)流場再次進(jìn)行數(shù)值模擬至收斂，計算內(nèi)壁面溫度Tw， in［i］。重復(fù)以上兩個步驟，直至滿足以下公式：

式（16）中，Tw， in［i］是當(dāng)前步第i次時推力室內(nèi)壁面溫度值，而Tw， in［i-1］是上一步第i-1次時推力室內(nèi)壁面溫度值，下標(biāo)“max”為發(fā)動機(jī)推力室計算域中相應(yīng)參數(shù)的最大值。

圖2　計算方法流程圖

3　計算結(jié)果與討論

本文采用不同的液滴參數(shù)對400N軌控發(fā)動機(jī)的工作過程進(jìn)行了數(shù)值模擬。

3.1Case 1數(shù)值模擬

3.1.1靜壓與靜溫

圖3、圖4顯示了當(dāng)前工況下的靜壓與靜溫，由圖可以看出，壓力在燃燒室內(nèi)最大，從噴管收斂段處開始減小，在噴管擴(kuò)張段繼續(xù)減小直至噴管出口。在拉瓦爾噴管中，在收斂段壓力開始降低直至噴管出口，燃燒室壓強(qiáng)為1.1MPa，比實驗所測（1.02MPa）略高，約為7.8%。發(fā)動機(jī)中心的溫度較高，發(fā)動機(jī)頭部的溫度相對較低，這是由于頭部有MMH液膜保護(hù)。

圖3　壓力分布

圖4　溫度分布

3.1.2MMH軌跡分布

圖5為液滴撞壁時采用反彈模型，即撞壁液滴全部反彈，且不考慮液膜時的MMH液滴軌跡分布，液滴全部反彈且全部蒸發(fā)，壁面無液膜；圖6為液滴撞壁時采用Bai-Gosman液膜撞壁模型且考慮液膜時，MMH液滴軌跡分布，MMH撞壁后一部分反彈，一部分形成液膜，但沒有飛濺，這是由于撞壁液滴的韋伯?dāng)?shù)不滿足Bai-Gosman液膜撞壁模型中液滴飛濺的條件。MMH撞壁形成的液膜可以對壁面起到熱保護(hù)作用。

圖5　反彈模型MMH軌跡

圖6　Bai-Gosman模型MMH軌跡

3.1.3壁面溫度分布

圖7為采用Stechaman半經(jīng)驗法獲得的壁面溫度分布，在當(dāng)前液滴參數(shù)下，計算的液膜長度小于實驗數(shù)據(jù)；發(fā)動機(jī)頭部溫度較低，未出現(xiàn)明顯的高溫區(qū)，這是由于MMH撞壁后在壁面形成了液膜，對頭部壁面起到了熱防護(hù)作用，使得溫度較低。

圖8為不考慮液膜冷卻時的壁面溫度分布，發(fā)動機(jī)頭部出現(xiàn)了局部高溫。這是由于MMH/NTO射流以相應(yīng)角度噴入推力室后，在壁面頭部存在MMH/NTO混合較均勻的位置，MMH/NTO發(fā)生化學(xué)反應(yīng)，而數(shù)值計算中未考慮液膜冷卻作用，導(dǎo)致該處壁面溫度較高。圖7、圖8對比分析表明，在發(fā)動機(jī)工作過程中，需考慮沉降型液膜冷卻，形成的液膜可對發(fā)動機(jī)壁面起到熱防護(hù)作用。

圖7　壁面溫度分布

圖8　壁面溫度分布

3.2增加MMH的質(zhì)量中徑（Case 2）

Case 2增加了MMH的質(zhì)量中徑。圖9為采用Stechaman半經(jīng)驗法獲得的壁面溫度分布，計算的液膜長度大于Case1，這可能是由于MMH質(zhì)量中徑增大，MMH蒸發(fā)變慢、撞壁概率增大，導(dǎo)致液膜長度增大；計算液膜長度稍大于實驗數(shù)據(jù)。圖10為液滴撞壁時采用Bai-Gosman液膜撞壁模型且考慮液膜時的MMH液滴軌跡分布情況，圖示MMH撞壁后一部分反彈，一部分形成液膜，這部分液膜可以對壁面起到熱防護(hù)作用。

圖9　壁面溫度分布

圖10　Bai-Gosman模型MMH軌跡

3.3增加NTO的質(zhì)量中徑（Case 3）

Case 3增加了NTO的質(zhì)量中徑。圖11為采用Stechaman半經(jīng)驗法時壁面溫度分布，在當(dāng)前液滴直徑下，其液膜長度比Case 2中低，這可能由于NTO質(zhì)量中徑增大，蒸發(fā)變慢，增加了NTO撞壁面的概率，增大NTO與壁面MMH液膜接觸的概率，接觸后，在附近發(fā)生化學(xué)反應(yīng)，消耗了形成液膜的MMH的流量，導(dǎo)致液膜長度比Case 2低，但計算的液膜長度與實驗中數(shù)據(jù)非常吻合。圖12為液滴撞壁面時采用Bai-Gosman液膜撞壁模型且考慮液膜時，MMH液滴軌跡分布，MMH撞壁后一部分反彈，一部分形成液膜，這部分液膜可以對壁面起到熱防護(hù)的作用。3.4液膜流率對液膜長度的影響

圖11　壁面溫度分布

圖12　Bai-Gosman模型MMH軌跡

Case 1、Case 2相比，僅有液滴的平均質(zhì)量中徑不同，其余參數(shù)均相同，而Case 3中液滴的最大直徑與Case 1和Case 2不同。采用Bai-Gosman液滴撞壁模型，采用自編程求解液滴撞壁后形成液膜時的流率并統(tǒng)計分析，分析結(jié)果如表2所示。

表2　液膜長度與液膜流率的分布

由表2可知，Case 1所得的液膜流率最小，Case 3所得的液膜流率其次，Case 2所得的液膜流率最大，Stechaman所得的液膜長度的順序與液膜流率成正比，即Case 2的液膜長度最長，Case 3的次之，Case 1的最小。這表明，液膜長度與液膜流量成正比，即液膜流量越大，液膜長度越大。

4　結(jié)束語

本文對小推力液體火箭發(fā)動機(jī)燃燒室內(nèi)的沉降型液膜冷卻進(jìn)行了數(shù)值模擬，得到以下主要結(jié)論：

（1）在Bai-Gosman液膜撞壁模型中，根據(jù)韋伯?dāng)?shù)將撞壁后液滴分為反彈、形成液膜和飛濺3種主要過程，可以較合理地數(shù)值模擬液滴撞壁后復(fù)雜狀態(tài)的變化過程。

（2）在沉降型液膜冷卻中，根據(jù)Bai-Gosman液膜撞壁模型，撞壁液滴只有在一定的韋伯?dāng)?shù)范圍內(nèi)，才形成液膜，而這部分流率是未知的，本文采用自編程求解液滴撞壁后形成液膜時的流率，并代入半經(jīng)驗公式方法中，液膜長度結(jié)果與實驗數(shù)據(jù)較吻合。

（3）本文在采用半經(jīng)驗公式方法分析液膜時，形成的液膜可以較好地對燃燒室的頭部起到熱防護(hù)作用，使頭部不出現(xiàn)明顯高溫。

（4）當(dāng)工況條件確定時，在沉降型液膜冷卻中，液膜長度主要由形成液膜的流率確定，液膜流率越大，液膜長度越大，對壁面的保護(hù)作用越好。

1Bai C X， Gosman A D. Development of methodology for spray impingement simulation ［J］. Sae Transactions， 1995，（104）: 21

2周紅玲， 楊成虎， 劉犇. 液體火箭發(fā)動機(jī)液膜冷卻研究綜述［J］. 載人航天， 2012， 18（4）: 8～13

3Ueda T， Tanaka H. Measurements of velocity， temperature and fluctuation distributions in liquid films［J］. International Journal of Multiphase Flow， 1975， 2（3）: 261～272

4Wasden F K， Dukler A E. Insights into the hydrodynamics of free falling wavy films［J］. AIChE Journal， 1989， 35（2）: 187～195

5Gater R A， L' Euyer M R， Warner C F. Liquid-film cooling， it's physical nature and theoretical analysis， AD627028［R］. Lafayette，Indiana: Jet Propulsion Center， Purdue University， 1965

6William M G. Liquid film cooling in rocket engines［Z］. 1991

7Stechman R C， Joelee O， Howel J C. Design criteria for film cooling for small liquid - propellant rocket engines［J］. Journal of Spacecraft and Rockets， 1969， 6（2）: 97～102

8Knab O， Preclik D， Estublier D. Flow field prediction within liquid film cooled combustion chambers of storable bipropellant rocket engines［R］. Aiaa Journal， 1998， 1998～3370

9Gaffney R L， White J A， Girimaji S S， et al. Modeling turbulent / chemistry interactions using assumed PDF method［Z］. Aiaa， 1992

10Seamans T F， Vanpee M， Agosta V D. Development of a fundamental model of hypergolic ignition in space-ambient engines ［J］. Aiaa Journal， 1967， 5（9）: 1616～1624

11林慶國， 楊成虎， 劉彝. 射流角度和壁面曲率對撞壁液膜的影響［J］. 國防科技大學(xué)學(xué)報， 2013， 35（2）: 17～21

1009-8119（2016）09（1）-0059-04