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

兗立文　許坤梅　王慧潔

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

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

兗立文許坤梅王慧潔

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

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

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

引言

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

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

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

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

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

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

1.2液滴運(yùn)動(dòng)控制方程

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

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

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

1.3Bai-Gosman液滴撞壁模型

Bai和Gosman［2］提出的模型將液滴撞壁過(guò)程詳細(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，θ是以弧度計(jì)算的液滴入射角。

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

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

二次液滴的速度方程為：

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

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

液滴的動(dòng)量方程為：

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

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

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

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

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

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

2　數(shù)值模擬

2.1計(jì)算參數(shù)

發(fā)動(dòng)機(jī)的噴注器采用雙組元離心式，向其內(nèi)噴嘴注入燃料甲肼（MMH），外噴嘴噴入氧化劑硝基-三唑-酮（NTO），氧燃推進(jìn)劑流率的混合比為1.667，全部流量為130g/s，初始溫度為300K。液滴的其余參數(shù)由前面計(jì)算得出。MMH的噴注速度Vinj，MMH為16.2m/s，NTO的噴注速度Vinj，NTO為2.6m/s，與對(duì)稱軸的夾角為55°。液滴進(jìn)入燃燒室時(shí)，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ù)詳見(jiàn)表1。

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

2.2計(jì)算網(wǎng)格

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

圖1　計(jì)算網(wǎng)格

2.3壁面邊界條件

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

2.4數(shù)值解法

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

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

圖2　計(jì)算方法流程圖

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

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

3.1Case 1數(shù)值模擬

3.1.1靜壓與靜溫

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

圖3　壓力分布

圖4　溫度分布

3.1.2MMH軌跡分布

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

圖5　反彈模型MMH軌跡

圖6　Bai-Gosman模型MMH軌跡

3.1.3壁面溫度分布

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

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

圖7　壁面溫度分布

圖8　壁面溫度分布

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

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

圖9　壁面溫度分布

圖10　Bai-Gosman模型MMH軌跡

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

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

圖11　壁面溫度分布

圖12　Bai-Gosman模型MMH軌跡

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

表2　液膜長(zhǎng)度與液膜流率的分布

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

4　結(jié)束語(yǔ)

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

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

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

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

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

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

2周紅玲， 楊成虎， 劉犇. 液體火箭發(fā)動(dòng)機(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林慶國(guó)， 楊成虎， 劉彝. 射流角度和壁面曲率對(duì)撞壁液膜的影響［J］. 國(guó)防科技大學(xué)學(xué)報(bào)， 2013， 35（2）: 17～21

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