低溫流體閃蒸噴霧研究

羅 盟，王鐵巖，張青松，王 楠，陳 慧，李文斌

(北京宇航系統工程研究所，北京 100076)

0 引言

液體在真空噴射過程中，由于壓力驟減(低于液體溫度對應的飽和蒸汽壓)，液體處于過熱狀態(熱力學亞穩定狀態)，在外界小擾動下，其迅速回歸到穩定狀態，并伴隨復雜的物理相變，產生“爆炸”式的霧化、汽化現象稱為閃蒸。對于閃蒸研究可以追溯到1960年，Brown等[1]開展了水和Freon-11介質的閃蒸試驗研究，并指出閃蒸存在臨界過熱條件。之后，Shepherd等[2]研究了過熱極限下液滴表面閃蒸的泰勒不穩定。Kita-mura等[3]、Park等[4]、Cleary等[5]研究提出了閃蒸的不同模式。Adachi等[6]在試驗基礎上提出了過熱蒸發速率的半經驗計算模型。Lamanna等[7]、Luo等[8]試驗研究了低溫流體的閃蒸特性，提出了判定完全閃蒸的無量綱參數。國內嚴俊杰等[9]、周致福等[10]、季璨等[11]、Guo等[12]在閃蒸換熱特性、噴霧場特性等方面開展了大量的工作，對閃蒸過程特性的認識具有重要的意義。

由于閃蒸過程涉及復雜相變過程，其理論并不成熟，目前仍缺少合適的數值計算模型。并且，在航天工程應用方面閃蒸現象具有較大的破壞作用，如歐洲阿里安火箭Aestus發動機曾出現過推進劑閃蒸“結冰”導致的點火燃燒問題[13]。鑒于此，本文在Zuo等[14]、Schmehl等[15]、Ramcke等[16]研究基礎上，建立了低溫流體閃蒸計算模型，開展了低溫流體的閃蒸過程氣-液兩相流仿真研究，并用試驗結果對其進行了驗證，為航天工程應用提供參考。

1 閃蒸噴霧數學物理模型

噴霧閃蒸可認為是由大量的單個液滴閃蒸組成。液滴閃蒸是一個涉及相變、傳熱傳質的復雜瞬態問題。閃蒸發生時，液滴的傳熱傳質機制主要包括內部過熱引起的核沸騰傳熱和外部的對流、熱傳導、輻射等傳熱。在較低的過熱條件下，外部傳熱具有重要作用，但在較高的過熱度下，內部過熱占據主導地位。因此，完整的閃蒸模型應同時考慮液滴內部過熱與外部傳熱。

1.1 過熱汽化換熱

研究推進劑噴射閃蒸過程，考慮液滴完整真實的傳熱傳質過程是十分復雜并且耗費資源的。因此，本模型中將液滴內部核沸騰作簡化處理，采用Adachi等[6]提出的半經驗參數模型。該模型假設過熱引起的換熱具有對流換熱的形式，并通過修正有效傳熱系數來考慮液滴內部過熱影響。其模型如下

(1)

(2)

式中，αf為內部有效傳熱系數，Ap為液滴表面積，Tp為液滴溫度，Tb為液滴飽和溫度，L(Tb)為飽和溫度下的汽化潛熱,ΔT=Tp-Tb。

1.2 外部導熱和對流換熱

在較低的過熱條件下，閃蒸過程中液滴外部氣體導熱和對流換熱占據重要的作用。在液滴蒸發過程準靜態、蒸發過程中始終保持理想球形，和液滴內部溫度均勻分布的假設下，Zuo等[14]給出了液滴導熱和對流換熱產生的蒸發質量流率

(3)

式中，r0為液滴初始半徑，λ為導熱系數，h為氣體焓值，cp為液滴比熱，Nu*為考慮沸騰的修正努塞爾數，其表達式為

(4)

Respe和Prspe為特定溫度Tspe下的雷諾數和普朗特數。BT為Spading數，其表達式為

(5)

氣-液混合熱力學參數和輸運特性參數均用特定溫度計算，該特定溫度采用Sparrow和Gregg提出的1/3規則[17]計算

(6)

式中，T∞為環境溫度，Tb為飽和溫度。

1.3 輻射換熱

因為低溫液體溫度較低，本文考慮了外界環境輻射對液滴表面蒸發的影響。由輻射換熱引起的液體蒸發質量流率為

(7)

式中，ε為液滴表面發射系數，本計算取0.96，σ為斯蒂芬-波爾茲曼常數，Ap為液滴表面積，Tamb為環境溫度。

1.4 液滴能量方程

閃蒸過程中，液滴內部溫度變化由過熱蒸發決定。液滴的能量方程為

(8)

1.5 液滴運動模型

運動顆粒軌跡的描述通過求解牛頓運動方程。在實際噴霧閃蒸霧化過程中，受到復雜的外力作用，如曳力、重力、壓力梯度引起的體積力、浮力等。本文忽略液滴之間的碰撞和融合，并認為氣動阻力是其受力主導因素。阻力模型采用修正的Schiller-Naumann模型[18]。液滴運動方程為

(9)

CD=

(10)

式中，vp和v分別為液滴和氣體速度矢量，BM為Spalding傳質數，CD為阻力系數，μ為流體的黏度，ρp為液滴密度，dp為液滴直徑，Re為雷諾數。

1.6 氣-液兩相相互作用

低溫推進劑真空噴射屬于兩相流動。本文在歐拉-拉格朗日框架下采用CFD-DPM方法對其進行仿真計算。通過求解三維非穩態可壓縮N-S方程來描述連續相(氣體)，通過求解牛頓運動方程追蹤離散相運動軌跡(液滴噴霧)，氣-液兩相之間的相互作用通過在連續相方程中添加相應的源項來考慮。三維非穩態N-S方程可參考文獻[19]。兩相之間的質量、動量、能量源項分別如式(11)～(13)

(11)

(12)

(13)

式中，Vcell為計算網格的體積。

1.7 物性參數

考慮到低溫推進劑閃蒸過程中溫度變化范圍較大(忽略壓力變化)，常物性參數計算帶來的影響較大，如圖1所示。因此本文基于NIST 9.1數據庫，采用多項式擬合來描述基于溫度變化的氣液物性參數(如密度、比焓，比熱、黏性、導熱率等)

圖1 模擬流體的物性Fig.1 Properties of the simulated fluids

φ(T)=a0+a1T+a2T2+a3T3+a4T4

(14)

通過以上分析建立了低溫流體閃蒸噴霧計算模型。本文將此模型植入FLUENT軟件中，當液滴溫度高于其周圍環境壓力對應的飽和溫度時，調用閃蒸模型；當溫度低于飽和溫度時，應用FLUENT的擴散控制蒸發模型(D2蒸發規律)。

2 試驗簡介

圖2給出了試驗裝置圖。試驗研究的是低溫推進劑經過單孔直噴嘴噴入低壓環境的過程。方形真空室為三面開玻璃窗結構，上端設計一個液氮儲罐，噴嘴、控制閥等浸泡在里面。試驗時，通過向小液氮罐加壓，控制其中液氮溫度。推進劑管路穿過小液氮罐，通過與液氮的換熱保證合適的推進劑噴射溫度。沿燃燒室中線軸線方向布置T型熱電偶(直徑1.0 mm)，每兩者間隔20 mm，用以測量噴霧沿軸線的溫度分布。燃燒室上布置壓力傳感器，用以測量真空室內壓力。噴嘴頭腔有熱電偶和壓力傳感器測口，用以測量噴射流體的溫度和壓力。噴霧形態采用高速紋影技術記錄。

圖2 試驗裝置圖Fig.2 Experimental setup

3 數值計算

3.1 計算方法

本研究基于CFD-DPM方法，在歐拉-拉格朗日框架下模擬氣-液(液滴)兩相流動。通過求解三維非穩態可壓縮雷諾平均N-S方程(U-RANS)計算連續相流場，通過求解牛頓運動方程追蹤離散相運動軌跡。氣-液兩相之間的相互作用是通過在連續相方程中添加相應的源項來考慮。采用兩方程標準k-ε湍流模型預測連續相湍流流動，采用隨機軌道追蹤模型來考慮連續相湍流效應對離散相的作用。仿真計算采用基于密度的隱式求解器，PISO壓力-速度耦合算法，動量方程、湍動能、湍流耗散率方程對流項均采用二階迎風格式離散，時間項采用一階隱式。

3.2 計算模型及邊界條件

計算模型及邊界條件如圖3所示。為節省計算資源，取模型1/4為計算域，采用對稱邊界條件。連續相采用氣相速度入口，壓力出口，絕熱無滑移壁面邊界條件。離散相壁面邊界采用彈性碰撞的反射邊界。入口離散相液滴大小以Rosin-Rammler分布規律給定。前期對計算空間離散和時間離散進行了敏感性分析。計算網格約180萬，連續相時間步長1×10-5s，離散相追蹤時間步長1×10-5s。

圖3 計算域、網格和邊界條件Fig.3 Computational domain, mesh and boundary conditions

4 結果與分析

本節通過液氮和液氧閃蒸的試驗結果驗證了仿真模型。由于液氧和液氮的閃蒸規律具有一定相似性，4.3～4.5節以液氧為例對仿真結果進行了詳細的分析。

4.1 噴霧形狀分析

圖4和圖5分別給出了仿真計算的液氧、液氮閃蒸噴霧中心截面上形態與試驗紋影對比圖。右半邊紋影圖中的黑色線條是熱電偶，直徑為1.0 mm。頂端噴嘴出口直徑為0.5 mm，可視化窗口寬度約為52 mm。試驗條件如表1所示。

表1 計算工況條件

由圖4～圖5可以看出，仿真計算較好地預測了液氧、液氮的閃蒸噴霧形態。該噴霧形態明顯不同于典型的射流噴霧(小錐形噴霧角)。閃蒸噴霧在噴嘴出口附近表現出很大的噴霧角和徑向穿透距離，這是由于低溫推進劑在噴入低壓環境的劇烈相變導致。

(a)仿真結果 (b)試驗紋影圖圖4 液氧噴霧Fig.4 LOx spray contour

(a)仿真結果 (b)試驗紋影圖圖5 液氮噴霧Fig.5 LN2 spray contour

以液氧工況為例，液氧噴入低壓環境后，在噴嘴出口，由于壓力的突然降低，液氧處于過熱狀態，在外界壓力擾動下釋放過熱回到穩定狀態。在回歸穩定狀態過程中是氣-液傳熱傳質的強烈非平衡過程，過熱液體中會產生大量初始氣泡核。初始氣泡核可以是均質結核也可以是異質結核。由分子運動理論可知，液體中分子團能量的不均勻性導致了流體密度起伏，局部低密度區形成初始的空化核心(均質結核核心)，而異質初始結核核心主要由固體接觸表面或者液體內部雜質提供。在過熱液體中形成一個半徑為r的氣泡核心后引起系統的自由焓的變化(等于形成該氣泡外界所需做出的功)為[20-21]

(15)

式中，σ為氣泡表面張力，Δμ為過熱液體與氣體的化學勢差，第一項恒為正，表示過熱液體生成氣泡需要做表面功。

隨著r值增大，系統的自由焓變隨氣泡大小變化曲線(紅線為均質結核情況，綠線為異質結核情況)如圖6所示。可以看出，該曲線存在最大值，其對應的氣泡半徑為臨界半徑。當氣泡半徑大于臨界半徑時，生成氣泡引起的系統自由焓變是減小的(生成氣泡所需的功減小)，這有利于氣泡的生長過程。在此情況下，大量氣泡生成，產生劇烈的核沸騰，氣泡的進一步生長使噴射流體破碎、霧化，伴隨著劇烈汽化，從而形成了如圖6所示的“爆炸式”噴霧形態。

圖6 自由焓與氣泡大小關系Fig.6 Free enthalpy versus sbubble size

4.2 噴霧溫度分析

閃蒸現象的一個重要研究對象是溫度分布特性。圖7和圖8分別給出了液氧、液氮噴霧軸線溫度分布仿真計算和試驗結果。圖中紅色虛線是噴霧軸線上離散相液滴溫度，紅色實線是軸線上連續相氣體的溫度分布。由于試驗中熱電偶測量的是熱電偶端頭局部區域內噴霧的溫度，即氣-液“混合團”的溫度。因此，本文對數值計算的液滴和氣體溫度采用1/3規則進行了平均，以此代替計算的局部氣-液“混合團”溫度，如圖7和圖8 中星號曲線所示?？梢钥闯觯抡嬗嬎憬Y果與試驗吻合較好，這也驗證了上述建立的閃蒸模型的正確性。由圖看出，低溫液體噴霧在噴嘴出口附近(Z<20 mm或Z/d<40)溫度由111 K急劇下降至85 K，降幅達26%，隨后在流場下游區域溫度分布平緩，降幅僅約5.8%。這間接預示了低溫液體在噴入真空燃燒室后，閃蒸現象隨即發生，并在噴嘴出口附近區域基本完成了過熱液體的閃蒸過程。

圖7 液氧噴霧溫度沿軸向分布Fig.7 LOx temperature distribution along the spray centerline

圖8 液氮噴霧溫度沿軸向分布Fig.8 LN2 temperature distribution along the spray centerline

4.3 噴霧粒徑變化分析

圖9和圖10分別給出了氧液滴粒徑D32(索泰爾直徑)沿噴霧軸向和徑向分布規律。圖9中的散點是在55 ms時刻軸線附近追蹤的顆粒大小分布。紅色曲線是0～55 ms時間段內的索泰爾直徑沿軸線分布的統計平均結果?？梢钥闯?，在噴嘴出口附近，氧液滴尺寸呈現出非常大的梯度分布規律，隨后液滴尺寸變化逐漸平緩。在0～20 mm范圍內，液滴尺寸從大約9.6 μm減小到8.5 μm，其減小幅度占整個計算域內減小幅值的70%。粒徑變化與溫度變化規律一致，這驗證了過熱閃蒸主要發生在噴嘴出口附近。圖中綠色點劃線是對散點圖進行了數據擬合，結果表明閃蒸過程中液滴大小較符合指數變化規律。若假設液滴在徑向各個截面的分布規律相似，則通過粒徑的軸向分布規律可以大致估算整個計算域內的蒸發量約占噴注量的34%。

圖9 液氧噴霧粒徑沿軸向分布Fig.9 LOx droplet size distribution along the spray centerline

圖10 液氧噴霧粒徑沿徑向分布Fig.10 LOx droplet size distribution in the spray radial direction

圖10給出離噴嘴出口不同軸向距離截面處，液滴大小沿徑向分布規律?？梢钥闯?，在靠近噴霧中心區域，液滴在各個不同截面的大小基本均勻。在同一截面上，噴霧外圍液滴尺寸沿徑向緩慢減小然后急劇增大，即較大的液滴聚集在噴霧外圍，這是由于大直徑液滴的慣性較大，受周圍連續相流場的影響較小，因此其更多地流向噴霧外圍。Munnannur等[22]以及Lasheras等[23]的兩相流噴霧試驗和計算也得到了類似的粒徑大小徑向分布規律。

4.4 噴霧脈動速度分析

在追蹤液滴過程中，液滴的運動由液滴和連續相氣體的相互作用決定。因此，液滴的脈動速度場也從一方面反映了整個噴霧流場的湍流程度。式(16)給出了追蹤液滴的均方根脈動速度(RSM)

(16)

圖11給出了沿噴霧軸線不同截面處，液氧噴霧脈動速度徑向分布曲線。可以看出，氧液滴在噴嘴附近(如Z=10 mm)表現出較大的速度脈動(～4%)，并且沿徑向分布波動較大。隨著軸向距離增加，速度脈動大幅減小，在噴嘴下游較遠處(如Z>100 mm)，脈動速度小于1%，并且在徑向外圍分布均勻。這是由于噴嘴出口附近的劇烈閃蒸引起的連續相流場的較大脈動，進而作用于液滴，引起液滴較大的速度脈動；而隨著軸向距離增加，閃蒸急劇減弱，整個噴霧流場趨于平穩。

圖11 液氧噴霧脈動速度沿徑向分布Fig.11 LOx velocity fluctuation along the spray radial direction

4.5 傳熱傳質分析

由以上的分析可知，低溫流體在噴入低壓環境后，在噴嘴出口附近區域便發生了劇烈的閃蒸現象，導致大的噴霧霧化角和徑向擴張，溫度急劇降低，液滴粒徑迅速減小；在噴嘴下游，隨著距離增加，以上參數變化平緩。閃蒸的這些特性是由其傳熱傳質決定的，本節對閃蒸過程中的傳熱傳質影響進行了分析。

圖12給出了液氧噴霧沿軸線方向的蒸發速率變化規律。由圖可知，在液氧噴入低壓環境后，強烈的閃蒸發生，產生的蒸發流量在10-9kg/s量級。之后，在較小的距離內(Z<20 mm)，閃蒸蒸發流量迅速降至10-11kg/s量級。在此區域，對流和導熱引起的蒸發流量較低，在10-12～10-11kg/s量級。因此，在靠近噴嘴出口附近，閃蒸占據絕對主導地位，高于其他換熱至少2個數量級。隨著噴霧流向噴嘴下游，由于液滴的過熱度大幅減小，閃蒸引起的蒸發流量快速減小，而對流和導熱逐漸占據主導作用。由于試驗前的預冷過程，真空燃燒室周圍初始環境溫度較低(約220 K)。因此，整個過程中的輻射換熱較弱。從計算結果可以看出，其產生的蒸發流量較低，基本維持在10-13kg/s量級。

圖12 液氧噴霧蒸發速率沿軸向分布Fig.12 LOx sprays evaporation mass flow rates along the spray centerline

4.6 閃蒸固態相變

值得一提的是，在非常大的過熱條件下，試驗中觀察到了液氮噴霧復雜的氣-液-固相變。固態氮更多地聚集在第一個熱電偶上(距離噴嘴出口20 mm)，并且靠近噴霧中心附近，如圖13所示。這也從一方面驗證了本文數值模擬結果的合理性，即過熱蒸發在靠近噴嘴區域占據主導地位，從而導致該區域換熱量非常大，產生了復雜的相變過程。在液氧噴射試驗過程中，由于試驗設備可提供的真空條件受限(高于40 mbar)，而氧氣三相點壓力為1.46 mbar，遠低于試驗真空條件，因此并未觀察到固態氧。相比而言，氮氣三相點壓力約為125 mbar，較容易觀察其凝固相變過程。

圖13 閃蒸過程中氮凝固現象Fig.13 Nitrogen solidification during the flashing

閃蒸的強烈換熱引起的固態相變會給火箭動力系統帶來潛在的危害，在實際工程設計中應予以重視。如發動機高空二次點火過程中，推進劑噴射閃蒸可能會產生固態推進劑塊，導致點火條件(點火處局部氧/燃比)偏離設計工況，帶來點火可靠性問題；或大量的固態推進劑團聚集在發動機壁面附近，推力室點火時可能產生爆燃，引起燃燒不穩定，甚至造成燃燒室結構件破壞。

5 結論

本文建立了低溫流體噴射閃蒸計算模型，在歐拉-拉格朗日框架下采用CFD-DPM方法開展了低溫液體真空噴射兩相流仿真計算，并進行了試驗驗證。結果表明該模型可以很好地預測閃蒸噴霧特性，包括閃蒸噴霧形態和溫度特性。研究結果表明：

1)低溫流體噴入低壓環境后，在噴嘴出口處便發生閃蒸現象，呈現出很大的噴霧角。

2)噴霧溫度、液滴大小在噴嘴附近呈現非常大的梯度變化。在出噴口較短的距離，噴霧溫度便降至接近環境壓力對應的飽和溫度，隨后溫度變化緩慢。

3)過熱閃蒸主要發生在噴嘴出口附近區域(x/d～40)，在此區域過熱閃蒸占據絕對主導地位(高出其他換熱2個數量級)；噴嘴下游，對流換熱與熱傳導影響占據主導地位；由于環境溫度較低，整個過程中輻射換熱量較小(低于其他換熱1～2個數量級)。

4)閃蒸帶來的強烈換熱可能導致噴嘴附近發生復雜的氣-液-固相變，會對發動機可靠點火帶來潛在風險。