徑向孔形狀對針栓式噴注器液膜下漏率的影響

王凱，雷凡培，楊岸龍，楊寶娥，周立新

1. 西安航天動力研究所 液體火箭發動機技術重點實驗室， 西安 710100

2. 中國船舶工業集團有限公司，北京 100044

近年來使用針栓式噴注器的針栓式發動機成為實現航天器垂直軟著陸的一大技術熱點和亮點，這是實現運載器重復使用和外行星探索的關鍵[1]，也成為實現低成本、可重復使用運載火箭回收的優選方案之一。

針栓式發動機經過五六十年的研究與發展，已取得了豐碩的工程成果，如著名的阿波羅登月艙下降發動機(LMDE)[2-3]、SpaceX公司的Merlin 1D系列發動機[4-5]以及中國的月球探測器CE-3用的7 500 N針栓式發動機。然而與工程研制的輝煌成就相比，目前針栓式噴注器工作特性的理論研究尚不足，許多機理尚不清楚，相關的研究主要集中于噴霧與燃燒特性研究，例如Ninish等[6]等開展了不同工況下測量霧化角和液滴粒徑的試驗研究及理論預測，Kim[7]和Sakaki[8]等分別開展了噴注單元的點火特性和燃燒特性試驗研究。在針栓式發動機性能提升的過程中，通常通過提高阻塞率來提高燃燒性能，然而往往面臨最突出的問題之一就是針栓頭燒蝕。在不對針栓頭采取額外熱防護措施的前提下，通過提高阻塞率提升性能似乎與針栓頭的熱防護問題是一對矛盾。工程上希望在不降低燃燒性能的前提下，通過采取一定的有效措施達到針栓頭的熱防護。

目前工程上使用較多的有效方法主要有2種：一種是被動熱防護，即針栓頭采用如鈮鎢合金等耐高溫的材料制成或者噴涂耐高溫抗氧化涂層，然而受材料極限溫度限制，單獨使用該方式多數情況下仍難以抵抗高熱流的沖刷燒蝕。另一種是主動熱防護，即通過在針栓頭上開孔，將中心一路的推進劑從頭部噴出實現主動冷卻。為了促使噴出的推進劑能較好霧化，常采用的主動冷卻孔形式有直流小孔[9]、自擊式噴嘴[10]、小離心式噴嘴[11]等。這種主動冷卻方式本質上通過改善流強與混合比分布的方法取得良好的熱防護效果。然而主動冷卻額外使用了一種推進劑，會降低總體燃燒效率，且難以適應大范圍變推要求。這是因為在變推力工況下，主動冷卻孔面積不可調，與變推力工況下推進劑主路噴注面積可調不匹配，導致低工況下主動冷卻孔路分配的流量占比顯著增大，對應的燃燒效率也明顯降低，且變比越大，低工況下這種效應越顯著。低工況下主動冷卻路流量占比可超過30%左右。Vasques和Haidn在研究LOX/LCH4針栓噴注器時用中心氧化劑路10%的流量主動冷卻，實現了固定結構工況下針栓頭的有效熱防護[9]。

經過上述分析可見，主動冷卻方式比較適用于固定推力或推力變比不大的情況，可以達到良好的針栓頭防護效果，而不太適用于大范圍變推力針栓噴注器，故主動冷卻方式不能完全滿足工程上的需求，即在不降低燃燒性能的前提下，達到針栓頭的熱防護；然而工程上希望在不降低燃燒效率的前提下，達到針栓頭的熱防護，即使在大范圍變推力的工況下，燃燒效率也不會顯著降低。

為了解決上述問題，本研究打算通過合理改變或組織噴霧過程來改善針栓頭附近的熱環境，從而滿足工程需求。LMDE未采用主動冷卻，表明依靠合理設計組織噴霧場的方案是可行的。通過合理改變或組織噴霧過程來改善針栓頭附近熱環境的想法其實源于針栓式噴注器噴霧場結構特征。針栓式噴注器是通過徑向液束與軸向環形液膜呈90°垂直撞擊，從而使推進劑進行霧化混合[12-13]，如圖1所示。針栓式噴注器的針栓頭伸入燃燒室，位于形成的大噴霧錐結構中心。燃燒狀態下針栓頭下游流場中存在的大回流區會攜帶大量熱燃氣回流至針栓頭附近。如果設計不合適，很容易造成針栓頭燒蝕。根據前期對針栓噴注器噴霧場的研究，發現軸向環形液膜與徑向液束撞擊后繞徑向液束形成分叉流動[14]，從而形成一定的液膜下漏流量。合理地設計并利用這部分下漏流量可達到與主動冷卻孔類似的效果，卻不用額外消耗推進劑，且在變工況過程中，下漏流量也會隨著主路推進劑一起調節變化，故具有較好的大范圍變推力流量匹配特性。合理地設計下漏流量，不僅可增加針栓頭附近的流強分布，將回流區推離針栓頭一定距離，減弱熱燃氣回流強度；同時還可改變針栓頭下游區域混合比分布，使其偏離當量混合比，從而降低該區域燃燒場的溫度。

圖1 針栓式噴注器原理圖

合理設計下漏流量的關鍵在于預估并控制其大小。如果下漏流量偏小，仍然容易引起針栓頭燒蝕；如果下漏流量偏大，過多的液膜下漏又會引起霧化混合效果變差，影響燃燒效率。因此，在設計針栓噴注器時，合理設計并準確預估下漏流量占原軸向液膜流量的比例(即下漏率)是研究如何提高針栓頭抗燒蝕能力的關鍵。

針對以上問題，本文以準確預估下漏率為切入點，以平面針栓多噴注單元為研究對象，基于前期針對徑向圓孔液束建立的考慮液膜液束變形和多噴注單元間相互影響的下漏率模型，首先通過類比分析，理論推導建立徑向矩形孔的膜束各自相對變形模型；再推導建立對應的不同高寬比矩形孔的下漏率模型；然后通過試驗及數值仿真對不同徑向孔形狀的下漏率模型進行驗證；最后通過對比分析為徑向孔形狀選擇給出合理建議，并給出模型中的常系數供工程設計預估使用，這對從設計初期就考慮針栓頭的熱防護問題具有重要的指導意義。

1 數值物理模型

1.1 數值方法

1.1.1 流動控制方程

對于噴霧計算過程來說，數值求解氣液2種流體均采用的是不可壓縮的、變密度的、帶有表面張力的Navier-Stokes方程[15]。具體的控制方程為

(1)

動量方程：

(2)

式中：u=(u,v,w)為流體速度；ρ=ρ(x,t)為流體密度；p=p(x,t)為流場中的壓力；τ為黏性應力；σ為表面張力系數；δ(x-x′)為Dirac函數，表示表面張力項集中在界面上；κ和n分別為界面的曲率和法向方向；S(t)表示氣液相界面。

另外，本文計算主要關注撞擊變形及流動過程，采用2種液體和氣體分別追蹤的分三相Volume of Fluid(VOF)方法。為了分別識別兩路液體的變形運動過程，定義了2種液相體積分數來進行2種推進劑的分別追蹤計算，分別為液相1體積分數α1和液相2體積分數α2，對應得到的流體密度和動力黏性系數為

ρ=α1ρl1+α2ρl2+(1-α1-α2)ρg

(3)

μ=α1μl1+α2μl2+(1-α1-α2)μg

(4)

式中：ρl1、ρl2和μl1、μl2分別為液相1和液相2的密度和黏度；ρg和μg分別為氣相的密度和黏度。

每一相液體對應的體積分數輸運方程為

(5)

1.1.2 氣液相界面捕捉及表面張力

計算中氣液相界面捕捉采用兩相流計算中最為常用的VOF方法，該方法具有較好的守恒特性[16]，且氣液相界面重構精度較高[17]。液體表面張力的計算采用常用的CSF(Continuum Surface Force)方法，該方法是通過在式(2)的動量方程中加入體積力源項來實現[15]。表面張力的計算式為

(6)

1.2 計算模型

在針栓式噴注器下漏率的計算模型中選取3個 完全相同的噴注單元作為研究對象，計算域模型如圖2所示。兩路液體均采用速度入口邊界，計算域兩側為對稱面邊界，計算中壁面均為無滑移邊界，其余面為壓力出口邊界，背壓設置為大氣環境。計算中選用兩路介質均為水，環境氣體為空氣，采用Fluent中標準的k-ε湍流模型和標準壁面函數[15]。計算域網格劃分采用880萬全結構化網格，對撞擊點及液膜流動區域進行加密處理，最小網格尺寸達到約30 μm，時間步長達到10-4ms量級，這樣可以快速收斂，進而更好地捕捉相界面和統計下漏流量。

圖2 多噴注單元計算域模型

2 多噴注單元試驗裝置及試驗測量系統

2.1 多噴注單元試驗裝置

冷態試驗的噴注器設計借用了文獻[8]中平面針栓式噴注單元的設計理念，采用噴嘴結構可更換的方案設計了平面針栓多噴注單元試驗裝置。平面針栓多噴注單元由一個可更換的液膜生成部分和一個可更換的多液束生成部分組成，試驗件如圖3所示。液膜生成部分可以產生一定厚度的平面液膜來模擬軸向液膜路，平面液膜足夠寬，設計中多孔噴注單元選取液膜寬度為20 mm；液束生成部分可以產生一定孔型型面和一定尺寸的射流。試驗中液膜和液束兩路單獨供應，通過分別改變液膜和液束兩路流量及更換不同結構尺寸的液膜和液束構件，從而形成不同的工況[18]。另外，原則上采用3個噴注孔即可模擬多噴注單元間的相互影響，然而試驗件中液束構件采用了9個徑向孔，這主要是為了增加展向寬度，便于下漏流量收集，同時也為了提高試驗中流量控制及下漏流量收集測量精度。

圖3 液液平面針栓多噴注單元試驗件結構

2.2 試驗測量系統

試驗系統簡圖如圖4所示，下漏流量采用在噴嘴下方采用容器收集的方案，試驗中使用工裝嚴格定位收集容器的位置，具體收集的噴霧場下漏區域與數值仿真統計的方式一致，即與軸向夾角為10°范圍內的液體流量，如圖5所示。在噴霧場達到穩定狀態時，通過收集一定時間的液體進行稱重來獲得收集的流量。試驗中分別對比了收集30、40、50、60、90 s的結果，分析發現收集時間超過50 s以上，換算得到的收集流量相差不多，相對誤差5%以內，收集30 s的結果與其有一定差別，故試驗中統一選取收集60 s進行試驗，且對于每個工況重復收集2～3次，取試驗平均值進行統計分析獲得下漏流量。

圖4 下漏流量收集試驗系統

圖5 試驗結果獲取下漏率示意圖

2.3 試驗數據處理

通過收集一定時間t內的液體，并稱重得到其質量為m，然而試驗中軸向液膜是具有一定寬度的，邊緣處存在無效撞擊液膜，直接無法測得有效的下漏流量。假設液膜在整個W=20 mm寬度內均勻，考慮到兩側的邊緣效應，統計時通過換算去除兩側的影響，如圖6所示，試驗件選用9孔也是為了降低邊緣效應的影響。于是有效下漏流量為8l+d寬度范圍內的流量，換算得到的下漏率計算公式為

圖6 下漏流量換算示意圖

(7)

3 理論建模分析

3.1 徑向圓形孔的下漏率基本定義及相關分析

對于針栓式噴注器，液膜與液束撞擊后部分液膜繞過液束會從相鄰液束孔之間下漏，形成液膜下漏流量。下漏過程如圖7所示，根據圖7中關系可以得到液膜下漏流量為

圖7 液膜/液束撞擊下漏過程示意圖

(8)

式中：ρ1為軸向路液體密度；u1為軸向路液膜速度；h為軸向路液膜厚度；w為與液束發生有效撞擊的液膜寬度。對應的阻塞率CBF=d/l。

(9)

式(9)是基于液束不變形建立的，實際中液束受撞擊后會發生變形，如圖7中虛線所示。隨著撞擊過程中液束的展向變形增大，更多的液膜流量被液束撞擊帶離軸向，實際液膜下漏流量及下漏率會較式(8)和式(9)預估的偏小，且液束變形越大，偏小程度越大。該影響主要體現在阻塞率上，由于液束的變形，展向寬度增加，實際阻塞率C′BF比幾何阻塞率CBF更大。

于是用實際阻塞率替換式(9)中的幾何阻塞率可得到考慮液束變形的下漏率預估模型：

(10)

式中：實際阻塞率C′BF=d′/l，d′為液束變形后的展向寬度(即長軸直徑)。可以發現，d′/d描述的是液膜/液束撞擊過程中液束的變形，w/d描述的是液膜/液束撞擊繞流過程中液膜的變形。

從式(10)可以看出下漏率與阻塞率CBF、w/d和d′/d有關。一旦確定了阻塞率CBF和液膜液束相對變形參數(d′/d和w/d)，便可預估出實際液膜下漏率。

3.2 徑向圓形孔的下漏率建模過程

(11)

從式(11)分析可以得出，w/d主要與動量比有關，即w/d=f(CMReff)，其在之前研究中已得到試驗換算結果及數值仿真結果的驗證，具有較好的準確性。另外，采用最小二乘法對試驗數據擬合得到了w/d與CMReff的關系式為

(12)

式中：a1為常系數。

通過數值仿真結果發現，液束受到液膜撞擊作用后從兩側開始變形，圓截面被壓扁，展向寬度逐漸增加，液束的撞擊前緣在液膜正壓和剪切作用下向下游移動，而后緣在液膜厚度范圍內幾乎未移動。因此，假設液束根部橫截面受到液膜撞擊作用后由圓形變成橢圓形，液束后緣位置保持不變，展向寬度d′增加(即為橢圓的長軸)，短軸可根據如圖8所示的幾何關系計算得到為d-h/tanθ，同時可假設液束流速大小保持不變，僅速度方向發生改變。

根據圖8的幾何關系及質量守恒可得

圖8 圓形孔液束根部變形過程模型示意圖

(13)

式中：θ′=(θ+π/2)/2，表示撞擊后液束中心速度方向為前后緣角度的平均值。

由式(13)可得液束的相對變形量為

(14)

實際的液束根部橫截面變形后雖然不是嚴格的橢圓形，與式(14)模型存在一定的差異，式(14)將變形后的橫截面當作橢圓處理，低估了變形量d′/d，但式(14)可以近似反映液束的相對變形d′/d與液膜/液束幾何特征參數h/d和有效動量比CMReff有關，可以近似地描述d′/d與h/d及CMReff的關系。式(14)也得到了如圖9所示的數值仿真結果的驗證，從添加流線的體積分數云圖上使用GetData軟件獲取邊界，測量獲取d′/d的定量數據。

圖9 添加流線的體積分數云圖

圖10 液膜撞擊液束根部受力變形分析

(15)

根據式(15)可對有效動量比很小時的情況做一定的修正。

以上述獲得的液膜相對變形量式(12)和液束相對變形量式(14)或(15)為基礎，將其代入下漏率表達式(10)中，并考慮常系數C0可得

(16)

另外，考慮到相鄰噴注單元之間的相互影響，引入相互影響系數來表征實際相互影響下的液束相對變形量與單個單元的比值，并假設其隨孔邊距(l-d)的增大而增大，取值范圍為0～1，即

(17)

相鄰多噴注單元之間的相互影響通過式(17)對液束的變形量進行修正，對液膜的變形量修正體現在指數a1中。修正后式(16)變為

Cleak=1-C0·CBF·

(18)

式(18)即為考慮多噴注單元間相互作用的液膜與徑向圓形孔液束撞擊的下漏率模型。

3.3 徑向矩形孔的下漏率建模過程

針對徑向孔為矩形的多噴注單元，類比上述徑向圓形孔，可得到類似的下漏率模型。液膜相對變形模型不變，只需將式(12)中的圓孔直徑d替換為矩形孔的寬度a，則式(12)變為

(19)

同樣，假設矩形孔受到液膜撞擊作用后橫截面仍近似為矩形，僅高寬比發生變化，即變扁，展向寬度a′增加，截面高度減小為b-h/tanθ，高寬比減小，如圖11所示。根據幾何關系計算可得到液束相對變形模型為

圖11 矩形孔液束根部變形過程模型示意圖

(20)

采用與徑向圓形孔同樣的思路，考慮小動量比下液束后緣移動效應的液束相對變形模型為

(21)

考慮多噴注單元間的相互影響后的下漏率模型變為

(22)

式中：阻塞率為CBF=a/l。

另外，液膜與徑向矩形孔液束相互作用過程中存在側邊效應的影響，參考文獻[6]中對矩形孔繞流側邊效應影響的表征，引入不同高寬比的矩形孔影響系數(1+CRb/a)，其中CR為對應的矩形孔系數，將其代入式(22)可得

(23)

式(23)即為液膜與矩形孔液束撞擊的下漏率模型。

4 模型驗證與結果分析

針對徑向圓形孔的液膜下漏率模型在之前的研究中已得到有效驗證[19]，分別對比了2種不同阻塞率結構下下漏率隨有效動量比的變化，均吻合較好。另外，在保證相同結構參數和有效動量比下，還對比了3種不同噴射速度下的下漏率，數值仿真、試驗結果與理論模型預測值均不隨噴射速度絕對值的變化而變化，進一步表明了下漏率理論模型的準確性。然而研究結果表明，徑向圓孔的下漏率隨著工況的降低會減小，這使得低工況下液膜下漏流量更低，造成低工況下熱防護風險增高，因此徑向圓孔似乎不太適合工況大范圍變化下的需求。

下面主要針對上述徑向矩形孔的下漏率模型開展不同有效動量比下的數值仿真算例驗證，研究對象為3種不同高寬比的矩形孔，孔的橫截面積均與d=1.0 mm的徑向圓孔相等。軸向液膜厚度為h=0.45 mm，徑向矩形液束孔寬度和高度分別取a=0.9 mm和b=0.9 mm、a=0.704 mm和b=1.15 mm、a=1.15 mm和b=0.704 mm，對應的高寬比b/a分別為1.0、1.63和0.61。相鄰孔中心間距均選取為l=2.0 mm，此時對應的阻塞率分別為CBF=0.45，CBF=0.35和CBF=0.575，具體的算例參數如表1所示，徑向圓孔對應的算例參數如表2所示。

表1 3種不同高寬比的矩形孔在不同有效動量比下的算例參數

表2 徑向圓孔在不同動量比下的工況參數表(CBF=0.5)

與徑向圓孔時一樣，根據數值仿真結果獲取下漏流量及下漏率的示意圖如圖12所示，提取與軸向夾角為10°范圍內的軸向液膜流量認為是下漏流量。為了減小采樣誤差，分別選擇6個橫截面位置作為采樣截面，即x0、x1、x2、x3、x4和x5，其中x0為徑向矩形孔液束中心所在截面位置，其余位置依次間隔1 mm。(注：此處選取10°夾角范圍與徑向圓形孔是一致的，這樣可以保證定量對比評估的準確性，這也是根據液膜下漏后的流場結構確定，且這樣選取可使得6個采樣截面的下漏流量接近，也與x0截面的下漏流量參考值相近)。試驗收集的噴霧場區域與數值仿真完全一致，詳見前文2.2節，通過式(7)統計分析得到對應的下漏率。

圖12 數值仿真結果獲取下漏率示意圖

獲取3種不同高寬比矩形孔的下漏率隨有效動量比的變化如圖13所示，分別對比了x0截面的下漏率、6個截面的平均下漏率、試驗測量下漏率及式(23)考慮噴注單元間相互影響模型的下漏率。式(23)模型中系數根據數值仿真結果擬合獲得，從圖13中可以發現數值仿真及試驗結果與理論模型預測值吻合較好，下漏率模型具有較高的準確性。同時可以發現3種不同高寬比情況下的下漏率均顯著小于幾何下漏率，即1-CBF。在所研究的工況下，下漏率均小于0.35。另外，下漏率隨有效動量比增大而增大的趨勢均較緩，特別是有效動量比大于2～3時。與如圖14所示的相同橫截面積的d=1.0 mm徑向圓孔相比，下漏率顯著減小，表明在孔心間距一定時，對于相同的徑向孔面積，矩形孔的下漏率小于圓孔的。從模型中的系數a1變小及C0(1+CRb/a)變大也可以看出，這正好說明矩形孔時的軸向液膜繞流效應相對較弱，液膜相對變形量隨有效動量比變化范圍較小，對應的實際阻塞率和實際下漏率隨有效動量比變化也較平緩。然而由于不同高寬比的矩形孔存在側邊效應，故高寬比b/a越大，幾何阻塞率越小，軸向液膜被帶到徑向的流量越小，使得下漏率越大。另外，通過徑向圓孔與矩形孔的對比，也可以推論得到矩形孔前緣兩側倒角可使得下漏率在原基礎上增大。

圖13 不同高寬比矩形孔的下漏率隨有效動量比的變化

圖14 徑向圓孔下漏率隨有效動量比的變化規律(CBF=0.5)

綜上所述，可以發現在徑向孔橫截面積及流量等工況參數完全一樣的情況下，徑向孔的形狀對下漏率有顯著的影響，矩形孔的下漏率低于圓形孔的；矩形孔的高寬比越大，下漏率越大。在實際應用中選擇矩形孔更有利于控制下漏率，并可通過改變高寬比來控制下漏率；同時在變工況調節過程中，只要保證h/b的值基本不變，無論改變噴注面積還是噴注壓降，下漏率均保持變化不大，下漏流量也會隨著主路推進劑一起調節變化，故具有較好的大范圍變推力流量匹配特性。

5 結 論

為了研究徑向孔形狀對針栓式噴注器液膜下漏率的影響并對其進行準確預估，考慮液膜液束各自變形和多噴注單元間相互影響推導建立了徑向圓形孔和不同高寬比矩形孔的實際下漏率模型，并通過試驗及數值仿真對模型進行了驗證分析。

1) 通過類比分析，基于徑向圓孔相對變形模型，推導建立了徑向矩形孔的相對變形模型及下漏率模型，并考慮了不同高寬比的矩形孔的繞流側邊效應，得到的理論預估結果與數值仿真及試驗結果吻合較好，表明針對矩形孔建立的相對變形模型及下漏率模型具有較好的準確性。

2) 針對矩形孔的研究表明，下漏率除了與幾何阻塞率、有效動量比及h/b有關外，還與高寬比有關。3種不同高寬比情況下的下漏率均顯著小于幾何下漏率，即1-CBF；高寬比b/a越大，下漏率越大。下漏率隨有效動量比增大而增大的趨勢均較平緩；且下漏率顯著低于相同橫截面積的徑向圓孔的。

3) 綜合分析徑向圓孔和3種不同高寬比矩形孔的結果發現，在工況參數完全相同的情況下，徑向孔的形狀對下漏率有顯著的影響。在實際中選擇矩形孔更有利于控制下漏率，并可通過改變高寬比控制下漏率；同時在變工況過程中，下漏流量也會隨著主路推進劑一起調節變化，保持下漏率變化不大，故具有較好的大范圍變推力流量匹配特性。