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

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

0 引言

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

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

1 閃蒸噴霧數(shù)學(xué)物理模型

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

1.1 過熱汽化換熱

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

(1)

(2)

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

1.2 外部導(dǎo)熱和對流換熱

在較低的過熱條件下，閃蒸過程中液滴外部氣體導(dǎo)熱和對流換熱占據(jù)重要的作用。在液滴蒸發(fā)過程準靜態(tài)、蒸發(fā)過程中始終保持理想球形，和液滴內(nèi)部溫度均勻分布的假設(shè)下，Zuo等[14]給出了液滴導(dǎo)熱和對流換熱產(chǎn)生的蒸發(fā)質(zhì)量流率

(3)

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

(4)

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

(5)

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

(6)

式中，T∞為環(huán)境溫度，Tb為飽和溫度。

1.3 輻射換熱

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

(7)

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

1.4 液滴能量方程

閃蒸過程中，液滴內(nèi)部溫度變化由過熱蒸發(fā)決定。液滴的能量方程為

(8)

1.5 液滴運動模型

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

(9)

CD=

(10)

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

1.6 氣-液兩相相互作用

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

(11)

(12)

(13)

式中，Vcell為計算網(wǎng)格的體積。

1.7 物性參數(shù)

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

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

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

(14)

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

2 試驗簡介

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

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

3 數(shù)值計算

3.1 計算方法

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

3.2 計算模型及邊界條件

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

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

4 結(jié)果與分析

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

4.1 噴霧形狀分析

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

表1 計算工況條件

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

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

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

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

(15)

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

隨著r值增大，系統(tǒng)的自由焓變隨氣泡大小變化曲線(紅線為均質(zhì)結(jié)核情況，綠線為異質(zhì)結(jié)核情況)如圖6所示?？梢钥闯?，該曲線存在最大值，其對應(yīng)的氣泡半徑為臨界半徑。當(dāng)氣泡半徑大于臨界半徑時，生成氣泡引起的系統(tǒng)自由焓變是減小的(生成氣泡所需的功減小)，這有利于氣泡的生長過程。在此情況下，大量氣泡生成，產(chǎn)生劇烈的核沸騰，氣泡的進一步生長使噴射流體破碎、霧化，伴隨著劇烈汽化，從而形成了如圖6所示的“爆炸式”噴霧形態(tài)。

圖6 自由焓與氣泡大小關(guān)系Fig.6 Free enthalpy versus sbubble size

4.2 噴霧溫度分析

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

圖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(索泰爾直徑)沿噴霧軸向和徑向分布規(guī)律。圖9中的散點是在55 ms時刻軸線附近追蹤的顆粒大小分布。紅色曲線是0～55 ms時間段內(nèi)的索泰爾直徑沿軸線分布的統(tǒng)計平均結(jié)果?？梢钥闯?，在噴嘴出口附近，氧液滴尺寸呈現(xiàn)出非常大的梯度分布規(guī)律，隨后液滴尺寸變化逐漸平緩。在0～20 mm范圍內(nèi)，液滴尺寸從大約9.6 μm減小到8.5 μm，其減小幅度占整個計算域內(nèi)減小幅值的70%。粒徑變化與溫度變化規(guī)律一致，這驗證了過熱閃蒸主要發(fā)生在噴嘴出口附近。圖中綠色點劃線是對散點圖進行了數(shù)據(jù)擬合，結(jié)果表明閃蒸過程中液滴大小較符合指數(shù)變化規(guī)律。若假設(shè)液滴在徑向各個截面的分布規(guī)律相似，則通過粒徑的軸向分布規(guī)律可以大致估算整個計算域內(nèi)的蒸發(fā)量約占噴注量的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給出離噴嘴出口不同軸向距離截面處，液滴大小沿徑向分布規(guī)律。可以看出，在靠近噴霧中心區(qū)域，液滴在各個不同截面的大小基本均勻。在同一截面上，噴霧外圍液滴尺寸沿徑向緩慢減小然后急劇增大，即較大的液滴聚集在噴霧外圍，這是由于大直徑液滴的慣性較大，受周圍連續(xù)相流場的影響較小，因此其更多地流向噴霧外圍。Munnannur等[22]以及Lasheras等[23]的兩相流噴霧試驗和計算也得到了類似的粒徑大小徑向分布規(guī)律。

4.4 噴霧脈動速度分析

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

(16)

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

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

4.5 傳熱傳質(zhì)分析

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

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

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

4.6 閃蒸固態(tài)相變

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

圖13 閃蒸過程中氮凝固現(xiàn)象Fig.13 Nitrogen solidification during the flashing

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

5 結(jié)論

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

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

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

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

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