小推力液體火箭發動機沉降型液膜冷卻液膜長度研究

兗立文　許坤梅　王慧潔

（北京航空航天大學宇航學院，北京　100191）

小推力液體火箭發動機沉降型液膜冷卻液膜長度研究

兗立文許坤梅王慧潔

（北京航空航天大學宇航學院，北京100191）

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

液體火箭發動機，液膜冷卻，液膜長度，液滴撞壁，數值模擬

引言

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

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

1　物理數學模型

1.1氣相湍流反應流的N-S方程

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

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

1.2液滴運動控制方程

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

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

其中，Re為液滴的雷諾數。

1.3Bai-Gosman液滴撞壁模型

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

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

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

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

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

借用氣象學中的相關數據擬合確定總的二次液滴數

二次液滴的速度方程為：

參考雨滴尺寸與速度的關系式推導出液滴的速度比為：

上式表明，若破碎后的液滴尺寸與速度呈負相關關系，尺寸越大，速度就越小。

液滴的動量方程為：

式中，Cf是摩擦系數，取值為0.6～0.8，θ1和θ2是相應的角度。

1.4Stechaman半經驗方法

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

Stechaman半經驗公式方法計算液膜長度為式（13）。

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

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

2　數值模擬

2.1計算參數

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

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

表1　噴霧噴注參數

2.2計算網格

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

圖1　計算網格

2.3壁面邊界條件

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

2.4數值解法

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

式（16）中，Tw， in［i］是當前步第i次時推力室內壁面溫度值，而Tw， in［i-1］是上一步第i-1次時推力室內壁面溫度值，下標“max”為發動機推力室計算域中相應參數的最大值。

圖2　計算方法流程圖

3　計算結果與討論

本文采用不同的液滴參數對400N軌控發動機的工作過程進行了數值模擬。

3.1Case 1數值模擬

3.1.1靜壓與靜溫

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

圖3　壓力分布

圖4　溫度分布

3.1.2MMH軌跡分布

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

圖5　反彈模型MMH軌跡

圖6　Bai-Gosman模型MMH軌跡

3.1.3壁面溫度分布

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

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

圖7　壁面溫度分布

圖8　壁面溫度分布

3.2增加MMH的質量中徑（Case 2）

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

圖9　壁面溫度分布

圖10　Bai-Gosman模型MMH軌跡

3.3增加NTO的質量中徑（Case 3）

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

圖11　壁面溫度分布

圖12　Bai-Gosman模型MMH軌跡

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

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

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

4　結束語

本文對小推力液體火箭發動機燃燒室內的沉降型液膜冷卻進行了數值模擬，得到以下主要結論：

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

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

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

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

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

2周紅玲， 楊成虎， 劉犇. 液體火箭發動機液膜冷卻研究綜述［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］. 國防科技大學學報， 2013， 35（2）: 17～21

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