胡文虎,孫皖,牛璐,彭婧,蔣卓君,潘良明

(1. 重慶大學低品位能源利用技術(shù)及系統(tǒng)教育部重點實驗室,400044,重慶;2. 重慶大學核工程與核技術(shù)系,重慶,400044; 3. 中國空氣動力研究與發(fā)展中心,621000,四川綿陽)

凝結(jié)作為一種常見的基礎物理現(xiàn)象廣泛存在于自然界以及工程領(lǐng)域,尤其在噴管、透平機、旋流分離器、低溫風洞等場景中,因工質(zhì)氣流極快,往往為跨聲速流動甚至是超聲速流動,更易導致凝結(jié)現(xiàn)象的發(fā)生。隨著高速飛行器的發(fā)展,對風洞全尺寸雷諾數(shù)試驗需求愈發(fā)迫切,常規(guī)風洞受限于自身種種原因無法繼續(xù)提升雷諾數(shù),通過液氮噴霧降低風洞內(nèi)部溫度并利用氣化后的氮氣作為試驗介質(zhì)被認為是現(xiàn)有提高試驗雷諾數(shù)方法中的最優(yōu)方法[1-2]。在溫度較低的情況下,氮氣由于膨脹降溫容易引發(fā)均質(zhì)凝結(jié)現(xiàn)象,此外,在極限運行總溫工況下生成的冰晶為低溫氮氣提供凝結(jié)核心,發(fā)生異質(zhì)成核,會使得凝結(jié)現(xiàn)象在相對較高的溫度下即可發(fā)生。這兩種不同成核形式的凝結(jié)現(xiàn)象最終均會影響氣動特性預測的準確性。

盡管在二十世紀就有了較為完善的凝結(jié)成核理論[3-5]和液滴生長理論[6-7],但是在過去很長的一段時間里,有關(guān)凝結(jié)流動的相關(guān)研究大多是基于濕蒸汽[8-12]和天然氣[13-15]進行的。受限于氮氣分子間勢能理論知識的匱乏,針對低溫氮氣開展相關(guān)研究存在不確定性和困難性,目前基于低溫氮氣凝結(jié)開展相關(guān)理論研究和實驗研究相對較少[16]。Faro、Wegener等[17-19]借助小型超聲速風洞進行了有關(guān)氮氣凝結(jié)的實驗研究;Willmarth、Bhabhe等[20-23]分別基于不同的測試手段研究了氮氣微尺度成核機理或過冷凝特性,但這些研究中凝結(jié)溫區(qū)均在氮氣三相點(63.15 K)以下。Sun等[24-26]使用不同的凝結(jié)成核模型和液滴生長模型在液氮級溫區(qū)對拉瓦爾噴管中的自發(fā)凝結(jié)進行了數(shù)值模擬研究,并借助構(gòu)建的模型完成了氮氣自發(fā)凝結(jié)現(xiàn)象對氣動性能試驗的影響分析,但未考慮風洞內(nèi)存在雜質(zhì)導致異質(zhì)成核凝結(jié)的影響。異質(zhì)成核方面,針對透平膨脹機以及超聲速旋流分離器背景下開展的研究比較多,陳紅梅、鞠鳳鳴等分別對一維噴管[8]和二維葉柵[9-10]內(nèi)濕蒸汽的異質(zhì)成核現(xiàn)象進行了數(shù)值研究,對比自發(fā)凝結(jié)闡明了存在雜質(zhì)對凝結(jié)過程的影響。邊江等[13-15]采用經(jīng)典成核理論和Gyarmathy液滴生長模型建立了天然氣均質(zhì)成核凝結(jié)以及異質(zhì)成核凝結(jié)模型,分析了不同變量對異質(zhì)成核的影響。

綜上,目前國內(nèi)外有關(guān)濕蒸汽以及天然氣膨脹凝結(jié)流動的研究成果斐然?；陲L洞中出現(xiàn)的凝結(jié)問題,不少學者也針對低溫氮氣開展了相關(guān)研究,但是目前少有液氮級溫區(qū)下關(guān)于雜質(zhì)組分存在時氮氣發(fā)生異質(zhì)成核的凝結(jié)流動特性的研究,突出表現(xiàn)在沒有一個統(tǒng)一的模型可以描述低溫氮氣不同成核形式下的膨脹凝結(jié)現(xiàn)象?？紤]濕蒸汽、天然氣膨脹凝結(jié)流動和低溫氮氣膨脹凝結(jié)流動過程具有一定程度的相似性,且相關(guān)理論體系比較完備,因此本文在濕蒸汽、天然氣膨脹凝結(jié)流動研究的基礎上建立了適用于不同成核形式的低溫氮氣膨脹凝結(jié)模型,在Fluent軟件中依托拉法爾噴管進行了數(shù)值模擬,分析了入口溫度以及異質(zhì)核心數(shù)和半徑對凝結(jié)流動的影響,為研究低溫風洞內(nèi)存在雜質(zhì)致使不同成核形式的凝結(jié)流動影響提供了理論依據(jù)。

1 數(shù)學模型

1.1 控制方程

基于Fluent軟件,構(gòu)建了低溫氮氣膨脹凝結(jié)兩相流動模型?？紤]氮氣在高速膨脹工況下形成的液滴尺度通常是微米量級,質(zhì)量非常小,形成的液滴隨著氣相流動在極短的時間內(nèi)就可以達到氣相的速度,因此可近似認為氣液兩相間無速度滑移。氣液兩相控制方程分別為

(1)

(2)

(3)

(4)

(5)

為了使方程組封閉還需添加液滴數(shù)守恒方程、濕度連續(xù)性方程以及液滴半徑液滴數(shù)和相關(guān)濕度之間的關(guān)系式

(6)

(7)

(8)

(9)

(10)

式中:下標hom、het分別表示均質(zhì)成核、異質(zhì)成核;J為成核率,m-3·s-1, 表示單位時間單位體積內(nèi)生成的液滴數(shù);Y為帶液量即液相質(zhì)量分數(shù);N為核心數(shù)密度,kg-1;rini為異質(zhì)核心半徑的初始值,m;若僅考慮均質(zhì)成核現(xiàn)象時將所有下標為hom、het的變量取值為0。

1.2 相變方程

(11)

(12)

式中:qc為凝結(jié)系數(shù),一般為1;ξ為非等溫修正系數(shù);σ為表面張力,N·m-1;ms為單個分子質(zhì)量,kg;kB為玻爾茲曼常數(shù),1.380 648 8×10-23J·K-1;Tg為氣相表面溫度,K;hfg為汽化潛熱,J·kg-1;Prg為氣體的普朗特數(shù);γ為比熱容比;Kn為克努森數(shù);αth為熱適應系數(shù),一般取1;Tl為液相溫度,K;r為液滴半徑,m。

質(zhì)量源項可以表示為

Sm=mhom+mhet

(13)

1.3 湍流方程和氣體狀態(tài)方程

RNGk-ε模型在Standardk-ε模型上進行了修正,可以更好地擬合彎曲壁面、流線以及高應變率的流動,但是對于強旋流精度不夠?？紤]本文數(shù)值模擬選用的噴管壁面以及流線的彎曲程度結(jié)合計算成本,本文選用RNGk-ε模型。低溫風洞通過噴射液氮降低風洞溫度,同時利用汽化后氮氣作為工作介質(zhì)進行試驗,其氣體狀態(tài)方程采用理想氣體狀態(tài)方程,液相密度、動力黏度、導熱率等相關(guān)物性參數(shù)通過NIST查詢、通過UDF模塊編程定義并加載到Fluent軟件中進行計算。

2 求解計算

2.1 計算域網(wǎng)格劃分及模型合理性驗證

本文依托拉法爾噴管進行數(shù)值模擬,噴管型線采用等膨脹率設計。采用如圖1(a)所示二維圖形拉伸而來的噴管,在軸向上具有對稱性,不同軸向位置截面上流動狀況無明顯差異,為提高計算效率,可在劃分網(wǎng)格時將模型簡化為二維。使用二維L型拓撲結(jié)構(gòu)對噴管進行網(wǎng)格劃分。噴管結(jié)構(gòu)參數(shù)、網(wǎng)格劃分后的噴管及邊界條件,如圖1所示。在進行數(shù)值模擬前,進行了網(wǎng)格無關(guān)性驗證。根據(jù)圖2最終選用網(wǎng)格數(shù)為45 260的網(wǎng)格進行了后續(xù)的模擬計算。

(a)噴管結(jié)構(gòu)參數(shù)

圖2 網(wǎng)格無關(guān)性驗證Fig.2 Mesh independence

為了驗證模型的有效性,模擬再現(xiàn)了計光華[27]進行的一組實驗工況。圖3對比了實驗結(jié)果和模擬結(jié)果靜壓沿著噴管流道中心線的變化情況,可見本文所建立的模型預測到了實驗中由于凝結(jié)相變釋放潛熱導致流道靜壓分布出現(xiàn)偏移的情況,由此可以在一定程度上驗證模型的有效性。

圖3 流道中心線上靜壓分布的數(shù)值結(jié)果與實驗結(jié)果的對比Fig.3 Comparison of static pressure distribution on the centerline of the flow channel between numerical simulation results and experimental results

2.2 數(shù)值計算邊界條件及工況

入口為壓力入口,入口液滴數(shù)以及液相體積分數(shù)均設置為0;出口為壓力出口,進口壓力設置為0.4 MPa,進出口壓力比為4;壁面設置為無滲流、無滑移、絕熱壁面。初場設置為標準初始化,設置初始表壓、速度、溫度、湍流等參數(shù)。

數(shù)值計算工況見表1,其中工況1～4用以探究入口溫度對均質(zhì)凝結(jié)現(xiàn)象的影響,工況5、6、10、12、13、14用以探究異質(zhì)核心數(shù)密度大小對凝結(jié)現(xiàn)象的影響;工況7～11用以探究異質(zhì)核心半徑大小對凝結(jié)現(xiàn)象的影響。其中異質(zhì)核心數(shù)及半徑參考均質(zhì)成核凝結(jié)現(xiàn)象中自發(fā)凝結(jié)核心的數(shù)值計算得到。

表1 數(shù)值計算工況

2.3 Fluent基本求解設定

求解器采用Ansys-Fluent,成核模型、液滴生長模型以及氮氣、氮液相關(guān)物性參數(shù)通過UDF定義,液滴數(shù)守恒方程、濕度連續(xù)性方程通過UDS定義,湍流模型采用RNGk-ε模型,收斂殘差目標設置為1×10-5。

3 均質(zhì)成核及異質(zhì)成核模擬結(jié)果

3.1 均質(zhì)成核模擬

從圖4均質(zhì)成核凝結(jié)現(xiàn)象發(fā)生后噴管內(nèi)液滴半徑沿著流道中心線的變化情況可知,液滴進入生長階段時半徑量級處于10-7m,始終處于μm量級以下,質(zhì)量非常小。由此可以說明,在建立模型前期假設忽略氣液兩相間的速度滑移是合理的。

圖4 凝結(jié)現(xiàn)象發(fā)生后液滴半徑分布Fig.4 Distribution of droplet radius after condensation phenomenon occurs

選用工況1～4的模擬結(jié)果,分析入口溫度對均質(zhì)成核凝結(jié)流動的影響。隨著非平衡凝結(jié)現(xiàn)象的發(fā)生,噴管中的流動會產(chǎn)生一定的影響。圖5(a)是靜壓沿著噴管流道中心線的變化對比情況。相比于無凝結(jié)現(xiàn)象發(fā)生的入口溫度為120 K的工況,93、95、97 K 3個工況靜壓變化具有相似性;發(fā)生膨脹凝結(jié)現(xiàn)象組的中心軸線靜壓分布有少許偏移。這是由于非平衡凝結(jié)過程考慮了相變中的熱力學非平衡過程,在發(fā)生凝結(jié)的位置,相變產(chǎn)生的潛熱加熱了氮氣流。流道中心線馬赫數(shù)變化如圖5(b)所示,由于潛熱的釋放,在凝結(jié)現(xiàn)象發(fā)生過后的噴管下游區(qū)域,馬赫數(shù)的增加也變得稍微緩慢。成核率變化如圖5 (c)所示,隨著入口溫度的降低,最大成核率出現(xiàn)位置前移,凝結(jié)現(xiàn)象發(fā)生的位置會逐步向喉部靠近,在有限長的噴管內(nèi),為液滴生長提供了更多的空間,液滴生長階段更長,因此出口帶液量隨著溫度的降低有所增加,如圖5(d)所示。

(a)靜壓分布

3.2 異質(zhì)成核模擬

圖6為噴管內(nèi)氮氣膨脹異質(zhì)成核凝結(jié)流動參數(shù)分布。氮氣進入噴管后隨著噴管的收縮,壓力、溫度逐漸降低,過冷度逐漸增加。與均質(zhì)成核凝結(jié)現(xiàn)象不同的是,當過冷度剛大于0時就發(fā)生了異質(zhì)成核凝結(jié)現(xiàn)象,可以觀察到異質(zhì)成核帶液量云圖中液體質(zhì)量分數(shù)逐漸增加。這是由于異質(zhì)核心的存在,氮氣分子不需要突破自由能障形成凝結(jié)核心,而是可以直接簇擁在異質(zhì)核心上發(fā)生凝結(jié)現(xiàn)象。隨著氮氣繼續(xù)在噴管內(nèi)膨脹,仍舊發(fā)生了均質(zhì)成核現(xiàn)象,具體情況與前文所分析的僅發(fā)生均質(zhì)成核凝結(jié)流動時的現(xiàn)象相同,不同的是由于異質(zhì)成核凝結(jié)現(xiàn)象也會伴隨相變釋放的潛熱,導致過冷度降低。從成核率云圖可以看出成核率有所降低,最大成核率從不添加異質(zhì)核心的1.18×1020m-3·s-1降低到了4.09×1019m-3·s-1,最大凝結(jié)核心數(shù)也從不添加異質(zhì)核心的2.887×1015kg-1減小到了1.433×1015kg-1,由此可見異質(zhì)核心的添加會對均質(zhì)成核現(xiàn)象產(chǎn)生一定的抑制作用。

噴管流道中心線上成核率以及凝結(jié)核心數(shù)的變化如圖7(a)、圖7(b)所示,可以看出,異質(zhì)核心數(shù)為1.0×1013kg-1時,最大成核率以及最大凝結(jié)核心數(shù)密度分別為1.177×1020m-3·s-1、2.883 9×1015kg-1,與無異質(zhì)核心添加時的模擬結(jié)果1.18×1020m-3·s-1、2.887×1015kg-1基本相同。噴管流道中心線上基于異質(zhì)成核產(chǎn)生的帶液量的變化如圖7(c)所示,可以看出,此時異質(zhì)成核帶液量也基本為0,可以認為基本無異質(zhì)凝結(jié)現(xiàn)象發(fā)生。隨著異質(zhì)核心數(shù)密度不斷增大,均質(zhì)成核現(xiàn)象被抑制的強度逐漸增大。異質(zhì)核心數(shù)密度增大到1.0×1016kg-1時,成核率以及凝結(jié)核心數(shù)基本沒有增長,噴管流道中心線上基于均質(zhì)成核產(chǎn)生的帶液量的變化如圖7(d)所示,此時均質(zhì)成核帶液量基本為0,可以認為均質(zhì)成核基本被完全抑制。

(a)成核率

當異質(zhì)核心半徑為1.0×10-6m時噴管流道中心線上異質(zhì)成核帶液量的變化如圖8(a)所示,可以看出,異質(zhì)成核帶液量基本上為0,可認為在該異質(zhì)核心半徑下噴管內(nèi)的凝結(jié)成核形式主要是均質(zhì)成核。隨著異質(zhì)核心半徑的不斷降低異質(zhì)成核現(xiàn)象愈發(fā)明顯,均質(zhì)成核現(xiàn)象在一定程度上被抑制,反映在帶液量變化(圖8 (a)、圖8(b))上是異質(zhì)成核帶液量逐漸增大,均質(zhì)成核帶液量逐漸降低。當異質(zhì)核心半徑達到1.0×10-8m以下時,可以發(fā)現(xiàn)隨著異質(zhì)核心半徑的降低,液相參數(shù)基本上沒有變化。考慮發(fā)生這種情況的主要原因是因為本次數(shù)值模擬的5組工況異質(zhì)核心數(shù)密度均為5×1014kg-1,不足以完全抑制均質(zhì)成核現(xiàn)象。可以預測,當異質(zhì)核心數(shù)充足時,隨著異質(zhì)核心半徑的減小,均質(zhì)成核現(xiàn)象會被完全抑制。

(a)異質(zhì)成核帶液量

3 結(jié) 論

本文在濕蒸汽高速膨脹凝結(jié)研究的基礎上建立了適用于不同成核形式的低溫氮氣跨聲速膨脹凝結(jié)模型,并在Fluent軟件中依托拉法爾噴管進行了數(shù)值模擬,計算結(jié)果表明:

(1)所建立的數(shù)學模型預測了由于噴管內(nèi)凝結(jié)現(xiàn)象釋放潛熱加熱內(nèi)部氣流使得壓力出現(xiàn)偏移造成下游馬赫數(shù)增速變緩的現(xiàn)象;此外隨著入口溫度的降低,凝結(jié)發(fā)生的位置會逐步向喉部靠近,出口液相參數(shù)會增大。

(2)相比于均質(zhì)成核,異質(zhì)核心的存在降低了成核的自由能障,使得氮氣凝結(jié)在0 K過冷度附近時即可發(fā)生;當異質(zhì)核心半徑量級處在10-7以下時,隨著異質(zhì)核心數(shù)的增多,異質(zhì)凝結(jié)逐漸成為帶液量增長的主導過程,異質(zhì)核心數(shù)密度達到1015量級以上時,發(fā)展為完全異質(zhì)成核。