周芮，程光平，張浩，任楓，王舜浩，張小斌

（1 浙江大學(xué)制冷與低溫研究所，浙江杭州310027； 2 上海宇航系統(tǒng)工程研究所，上海201109）

引 言

某新型運(yùn)載火箭擬采用冷氦增壓技術(shù)，在國際上中國是第二個(gè)掌握這項(xiàng)技術(shù)的國家[1-2]。冷氦增壓技術(shù)應(yīng)用于使用低溫推進(jìn)劑一類的火箭中，氦氣瓶放置在液氧箱中，低溫氦氣進(jìn)入煤油箱增壓。利用高壓低溫下，氦氣密度大的特點(diǎn)，一個(gè)冷氦氣瓶置于液氧箱中的貯氣量約為同壓力、同容積下常溫氣瓶貯氣量的3倍。因此該技術(shù)可以大大減輕貯箱系統(tǒng)質(zhì)量，提高火箭運(yùn)載能力。根據(jù)新型運(yùn)載火箭設(shè)計(jì)方案，擬采用低溫氦氣在液面以下以直接鼓泡方式噴出，然后與煤油換熱升溫后進(jìn)入氣相空間，提高氣體溫度和增壓能力。但這種煤油貯箱增壓方法存在與冷氦氣流接觸的局部煤油溫度降低結(jié)冰的風(fēng)險(xiǎn)，結(jié)冰沉底將堵塞發(fā)動(dòng)機(jī)濾網(wǎng)，從而引起發(fā)動(dòng)機(jī)故障。另外煤油排出過程可能夾雜氦氣泡進(jìn)入發(fā)動(dòng)機(jī)，誘導(dǎo)發(fā)動(dòng)機(jī)泵發(fā)生氣蝕，使其不能正常工作。而且氣泡在煤油中上升速度慢，可能存在增壓延遲。針對這些可能的問題，本文對冷氦氣在貯箱中心單孔噴入和環(huán)向多孔噴入兩種工況（兩種工況冷氦氣流量相同）開展了計(jì)算流體力學(xué)（CFD）數(shù)值分析，分別建立了軸對稱和局部三維幾何模型。為驗(yàn)證數(shù)值模型，計(jì)算結(jié)果首先與基于原理性實(shí)驗(yàn)貯箱，冷氦氣在常溫水中的鼓泡實(shí)驗(yàn)觀察進(jìn)行了對比，并開展了煤油的數(shù)值計(jì)算和分析。

國內(nèi)研究者對于冷氦增壓系統(tǒng)展開了許多模擬和實(shí)驗(yàn)研究[3-5]。針對系統(tǒng)的研究，曾源華[6]分析了模擬火箭飛行過程中加熱20.4 K低溫氦的實(shí)驗(yàn)結(jié)果。邢力超等[7-8]建立了液氫溫區(qū)冷氦增壓系統(tǒng)，采用深低溫制冷機(jī)配合高壓低溫?fù)Q熱貯罐獲得20 K、21 MPa 冷氦氣源，同時(shí)使用管殼式加熱器模擬火箭上的換熱器性能，實(shí)驗(yàn)系統(tǒng)安全可靠性得到了證實(shí)。針對具體的增壓過程，白曉瑞[9]采用AMESim 軟件對液氧增壓系統(tǒng)進(jìn)行方案設(shè)計(jì)，并且比較了氦氣、氮?dú)夂涂諝獾脑鰤盒Ч?，最終結(jié)果表明氦氣的增壓效率最高。王承等[10]利用仿真軟件Sinda/Fluint基于集總參數(shù)法對冷氦氣瓶的加注過程和系統(tǒng)增壓過程進(jìn)行建模和計(jì)算分析，結(jié)果對氣枕壓力和氣瓶熱力狀態(tài)之間的關(guān)系進(jìn)行了深入探討。針對增壓系統(tǒng)中的換熱器研究，杜正剛等[11]對各工況下?lián)Q熱器的換熱性能和壓力損失展開了實(shí)驗(yàn)和模擬研究，結(jié)果表明冷氦氣流量增大會(huì)導(dǎo)致?lián)Q熱器出口溫度降低和壓力損失。趙耀中等[12]對比了三種換熱器，結(jié)果表明采用低溫氦氣-常溫氮?dú)獾膿Q熱器最好，原理恰當(dāng)，容易實(shí)現(xiàn)且便于操作。

計(jì)算流體力學(xué)（CFD）方法也被廣泛應(yīng)用于氦氣增壓低溫推進(jìn)劑貯箱（液氧和液氫）熱力學(xué)過程的研 究[13-19]。NASA Glenn 研 究 中 心 的Kassemi 等[20]基于Fluent 軟件開展了二維軸對稱液氫貯箱自增壓、循環(huán)去熱分層過程的CFD 數(shù)值研究。結(jié)果發(fā)現(xiàn)Lee模型的經(jīng)驗(yàn)系數(shù)對計(jì)算結(jié)果影響較小。即使Reynolds 數(shù)顯示流動(dòng)是完全湍流范圍，層流模型計(jì)算得到的壓力動(dòng)態(tài)變化比基于湍流模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果更一致,原因是湍流有效熱導(dǎo)過大。國外有研究者基于Fluent軟件模擬內(nèi)壁帶方形和半圓形肋片的圓柱貯箱的自增壓過程，目的是檢驗(yàn)肋片結(jié)構(gòu)對貯箱壓力動(dòng)態(tài)變化的影響，研究了方形肋片相對間隔以及肋片結(jié)構(gòu)的影響[21]。分析基于2D 軸對稱模型，兩相流基于VOF 模型，湍流采用二方程standard k-ε 模型。結(jié)果顯示半圓形肋片對減小壓力增加速率效果比方形肋片更好。國內(nèi)低溫推進(jìn)劑國家重點(diǎn)實(shí)驗(yàn)室聯(lián)合西安交通大學(xué)基于CFD 方法開展了液氧貯箱在加壓加注階段（頂部噴入高溫氧氣）以及發(fā)射過程最后級自增壓的研究[22-23]。所采用的兩相數(shù)值模型同樣基于2D 軸對稱及VOF 模型，湍流采用二方程standard k-ε 模型。系統(tǒng)研究了貯箱絕熱結(jié)構(gòu)外環(huán)境溫度對自增壓的影響，考慮了空間輻射的影響。研究結(jié)果顯示，環(huán)境溫度變化（253 ～313 K）對液氫增壓幾乎沒有影響。Zhou 等[24]用Fluent 基于VOF 方法，采用層流模型模擬了地面液氫儲罐的自增壓過程，探討了儲罐臨界壓力對蒸發(fā)過程的影響?？紤]到儲罐實(shí)際的結(jié)構(gòu)，模擬過程中理論推導(dǎo)了兩種儲罐表面（包裹泡沫和不包裹泡沫）是否結(jié)冰的熱力學(xué)過程，并通過UDF 實(shí)現(xiàn)并應(yīng)用于邊界條件。通過與實(shí)驗(yàn)數(shù)據(jù)的對比，驗(yàn)證了該模型的準(zhǔn)確性，深入探討了液氫的三個(gè)蒸發(fā)階段并且為低溫儲罐的設(shè)計(jì)提供了參考。

以上研究內(nèi)容都是對室溫氦氣進(jìn)入低溫貯箱氣相空間的增壓過程進(jìn)行數(shù)值研究，低溫氦氣進(jìn)入液相鼓泡增壓的數(shù)值研究報(bào)道較少[25]。對比兩種物理過程，前者氣液間滑移速度(slip velocity)可忽略不計(jì)，兩相流動(dòng)計(jì)算都是基于體積含量模型（VOF）。該模型為“單流體”模型，即假設(shè)氣液兩相流為一種混合流體，同時(shí)不考慮相間的相對速度。而后者低溫氦氣高速噴射進(jìn)液體，氣液間相對速度不能忽略，因此有必要構(gòu)建Euler-Euler 模型。相比于VOF模型，Euler-Euler 模型對兩相分別列出控制方程，并考慮相間相對速度，但同時(shí)也更復(fù)雜，對計(jì)算資源要求更高。本文基于Euler-Euler 模型對煤油貯箱中冷氦鼓泡增壓進(jìn)行了數(shù)值模擬，模擬結(jié)果與水箱的冷氦鼓泡實(shí)驗(yàn)進(jìn)行了對比。所得結(jié)果希望為某新型運(yùn)載火箭的煤油增壓方案設(shè)計(jì)提供參考。

1 數(shù)值模型

建模的貯箱直徑1.2 m，圓柱段高度1.8 m，兩端為橢球形封頭，直徑比1.6，底部有排出口。圖1 和圖2 分別給出了氦氣單孔鼓泡和環(huán)形鼓泡增壓CFD幾何模型及網(wǎng)格，其中單孔鼓泡模型的網(wǎng)格數(shù)為49386 個(gè)，中心孔半徑為15 mm，在貯箱內(nèi)部的氦氣輸送管由于管徑相對較小，沒有被建模。環(huán)形鼓泡模型的網(wǎng)格數(shù)為1428008 個(gè)，環(huán)向扇形角為30°，每隔14.4°有一個(gè)直徑3 mm 的小孔。在兩種鼓泡方式中，氦氣孔到封頭底部之間的距離均為0.575 m。

考慮到低溫氦氣與煤油之間的溫差和速度差，有必要構(gòu)建更全面復(fù)雜的Euler-Euler 兩相模型。同時(shí)為減小模型復(fù)雜性，本文通過煤油的溫度以及是否有時(shí)間結(jié)冰兩步來判斷結(jié)冰可能性，因此模型中不包括液固相變，相應(yīng)的控制方程如下。

質(zhì)量守恒方程

動(dòng)量守恒方程

式中，k 指氣相（G）或液相（L）。τD,k為氣液相間拖曳力，拖曳力的表達(dá)式為

能量守恒方程

圖1 單孔鼓泡增壓方案儲罐網(wǎng)格Fig.1 Mesh scheme of central injection from a core

圖2 環(huán)形鼓泡擴(kuò)散增壓網(wǎng)格Fig.2 Mesh scheme of circumferential injection from three cores

考慮到流場中既有高Re區(qū)域（氦氣噴出口），也存在低Re 區(qū)域（遠(yuǎn)離噴出口的區(qū)域），因此選擇SST k-ω 湍流模型。該模型綜合了k-ω 模型在近壁區(qū)計(jì)算的所有優(yōu)點(diǎn)，并在湍流黏度定義中考慮了湍流剪切應(yīng)力輸運(yùn)過程，適用更廣，可以用于帶逆壓梯度的流動(dòng)計(jì)算和翼型計(jì)算等。

計(jì)算中所有的物性數(shù)據(jù)都來自Refprop 9.1。對于氦氣：密度由理想氣體公式計(jì)算。壓力和溫度的影響很小，因此為常數(shù)cp=5193 J/(kg·K)。熱導(dǎo)率和動(dòng)力黏度為溫度的函數(shù)。本文首先建模了低溫氦氣在水中的鼓泡增壓過程，以與實(shí)驗(yàn)對比。計(jì)算中，水的物性都為常數(shù)，參考溫度為20 ℃。

煤油物性也都為常數(shù)，其中密度為805 kg/m3，比熱容為1967.8 J/(kg·K)，熱導(dǎo)率為0.149 W/(m·K)，動(dòng)力黏度為0.001 kg/(m·s)。

模擬中，單孔和環(huán)形鼓泡中氦氣和煤油總排放量分別為0.04 kg/s 和55.4 kg/s，其中環(huán)形鼓泡方案中的每個(gè)小孔流量為0.0016 kg/s。整個(gè)流體域初始溫度為280 K，絕對壓力為0.32 MPa。

基于Fluent18.1 進(jìn)行一階隱式非穩(wěn)態(tài)數(shù)值計(jì)算。迭代方案：壓力-速度耦合方法Coupled，相比于壓力和速度Separated 方案，Coupled 方案同時(shí)求解動(dòng)量方程和質(zhì)量守恒方程，同步收斂壓力和速度，因此在計(jì)算非穩(wěn)態(tài)時(shí)能夠更快地收斂。其他方程的離散格式設(shè)置如下：壓力求解，PRESTO!；動(dòng)量方程，二階迎風(fēng)格式，體積含量，一階迎風(fēng)格式；能量方程，一階迎風(fēng)格式。方程收斂標(biāo)準(zhǔn)：質(zhì)量守恒方程殘差RRR＜10-3，所有速度量殘差RRR＜10-3，能量殘差RRR＜10-6。

2 模型驗(yàn)證

圖3 冷氦環(huán)形鼓泡實(shí)驗(yàn)貯箱內(nèi)部Fig.3 Interior of ringlike holes scheme tank

圖4 冷氦在水中環(huán)形鼓泡實(shí)驗(yàn)快照Fig.4 Experiment moment during cold helium pressurization process in water tank

為驗(yàn)證方案可行性及數(shù)值模型，上海宇航系統(tǒng)工程研究所和上海交通大學(xué)黃永華課題組合作開展了90 K 氦氣在水中環(huán)形鼓泡實(shí)驗(yàn)，圖3 給出了內(nèi)部相機(jī)從上往下拍攝的實(shí)驗(yàn)裝置內(nèi)部，包含冷氦進(jìn)入管、熱電偶支架以及環(huán)形排氣管等，進(jìn)入管和環(huán)形管通過三根水平方向的旁管聯(lián)通。圖4給出了另一個(gè)攝像機(jī)從水平方向拍攝的環(huán)形排氣管附近冷氦氣剛開始排放時(shí)的現(xiàn)象。從圖4 可以看出，環(huán)形管上各個(gè)小孔流量并不均勻，靠近旁通管的小孔首先噴出氦氣，而后其他小孔噴出氣體。氦氣噴出小孔后，馬上與相鄰小孔噴出的氣體混合，形成大氣團(tuán)，攜帶液體上升，即存在明顯的氣液混合和攜帶現(xiàn)象。

圖5 給出了以小孔為中心，左右偏轉(zhuǎn)7.5°三個(gè)截面上典型的冷氦氣鼓泡氣液相分布，7.5°的偏轉(zhuǎn)截面處于兩個(gè)小孔的中間。由于三個(gè)截面都出現(xiàn)了氦氣，可見冷氦氣噴入水中后擴(kuò)散并與相鄰小孔的氦氣進(jìn)行了混合，與實(shí)驗(yàn)現(xiàn)象一致。另外布置的熱電偶溫度傳感器由于離環(huán)形管較遠(yuǎn)，實(shí)驗(yàn)過程沒有觀察到溫度明顯下降，與數(shù)值結(jié)果一致。實(shí)驗(yàn)過程所有管路上沒有觀察到結(jié)冰現(xiàn)象。

3 計(jì)算結(jié)果分析

3.1 單孔鼓泡

圖5 環(huán)形鼓泡狀態(tài)下模擬的不同截面上冷氦氣從水中鼓泡噴出結(jié)果（“1”指的是全氣相）Fig.5 Simulation results of cold helium pressurization process in water tank on three different sections (legend denotes volume fraction of gas)

圖6 單孔鼓泡時(shí)儲罐中的壓力隨時(shí)間的變化曲線Fig.6 History of ullage pressure in one hole scheme tank

圖6給出了中心軸上單孔鼓泡時(shí)貯箱氣相空間壓力隨時(shí)間的變化。考慮到實(shí)驗(yàn)中以水代替工質(zhì)，因此分別模擬了液相為水和煤油兩種工況，以定量說明水是否可以代替煤油得到一致的實(shí)驗(yàn)結(jié)果。首先，由于排出液體體積流量和氦氣進(jìn)入流量完全相同，因此兩者的壓力變化趨勢幾乎一致，但數(shù)值上水儲罐的壓力始終比煤油儲罐壓力稍小，原因可能是水罐氣相空間中氦氣溫度比煤油儲罐低。圖7給出了t=36 s時(shí)刻水儲罐和煤油儲罐氣枕空間的溫度云圖分布，從圖中可以明顯看到水儲罐氣枕空間內(nèi)，氦氣的平均溫度更低。冷氦氣從噴口出來后，與高溫液體進(jìn)行強(qiáng)烈混合對流和換熱。由于煤油密度比水?。?.8∶1），氦氣在煤油中的浮力比在水中更小，到達(dá)液面的時(shí)間更長，計(jì)算顯示，當(dāng)工質(zhì)為水時(shí)，最大速度為5.52 m/s，煤油中最大速度為5.24 m/s，因此煤油和氦氣混合和換熱更充分，導(dǎo)致氦氣到達(dá)液面后溫度更高。其他結(jié)果，如溫度、相分布等定性一致，因此下面只對氦氣/煤油數(shù)值結(jié)果進(jìn)行分析。

圖8 給出了氦氣/煤油單孔鼓泡在t=4,16,32 s時(shí)刻的溫度分布云圖。 由圖可見，由于(ρcp)氦氣?(ρcp)煤油，氦氣在與液體換熱后溫度很快就上升到了高溫，自始至終低溫區(qū)域被限制于氦氣入口處。氦氣相對于煤油高速噴出，因此是強(qiáng)制對流換熱，同時(shí)存在強(qiáng)烈的湍流動(dòng)量交換，在出口附近被冷卻的煤油被氦氣帶走，同時(shí)更遠(yuǎn)的熱流體流過來補(bǔ)充，溫度場達(dá)到動(dòng)態(tài)平衡。此過程中，由于熱傳導(dǎo)引起的熱擴(kuò)散相對于對流換熱影響較小，因此，溫度場沒有隨著時(shí)間延長擴(kuò)散。同時(shí)，受煤油側(cè)的流動(dòng)和壓力場擾動(dòng)，以及強(qiáng)烈湍流的影響，溫度場分布并不能保持不變的分布。

圖7 t=36 s時(shí)刻水儲罐和煤油儲罐氣枕空間的溫度云圖分布Fig.7 Temperature contours of gas phase inside water/kerosene tank at t=36 s

圖8 氦氣/煤油單孔鼓泡在t=4,16,32 s時(shí)刻的溫度分布云圖Fig.8 Temperature contour of helium/kerosene in central injection scheme tank at t=4,16,32 s

圖9 在t=4 s時(shí)刻x=0.664 m處的溫度和液相含量分布曲線Fig.9 Temperature and liquid fraction lines at x=0.664 m and t=4 s

圖9給出了在t=4 s時(shí)刻在噴口上方0.006 m（即x=0.664 m 處）沿水平方向（如圖8 實(shí)線所示）的溫度和相含量分布。要使結(jié)冰現(xiàn)象發(fā)生需要滿足三個(gè)必要條件：①有液體存在；②溫度低于固化溫度（這里煤油固化溫度為T=213.15 K）；③時(shí)間足夠長。從圖中可以看出，滿足①和②兩個(gè)條件的區(qū)域在0.015 m≤y≤0.019 m，寬度約4 mm，該區(qū)域?qū)?yīng)的液相含量約為0～0.04。因此該區(qū)間內(nèi)的液滴是否結(jié)冰決定于第③個(gè)條件。當(dāng)噴口附近的煤油液滴溫度降到固化溫度213.15 K（-60℃），要使其進(jìn)一步固化，則溫度需要從液滴外表面開始到中心進(jìn)一步降低（考慮到一定的過冷度）。忽略液滴內(nèi)部對流換熱，從外表面到球心為一維導(dǎo)熱，則根據(jù)傅里葉系數(shù)（Fo）

3.2 環(huán)形鼓泡

圖10 給出了環(huán)形鼓泡時(shí)貯箱氣枕壓力隨時(shí)間的變化，為對比，圖中同時(shí)給出了中心孔鼓泡時(shí)的氣枕壓力變化。由圖可見，相比于中心孔進(jìn)氣，環(huán)形進(jìn)氣時(shí)氣枕壓力存在明顯的滯后，且隨著時(shí)間增長壓差越來越大。原因是環(huán)形進(jìn)氣時(shí)每個(gè)小孔流速更小，因此到達(dá)氣枕的時(shí)間更長，同時(shí)存在與液體的混合與旋渦（圖5 和圖11），部分氣體長時(shí)間滯留在液相中，所以環(huán)形進(jìn)氣時(shí)氣枕壓力上升緩慢。若要讓儲罐達(dá)到相同的壓力，在采用環(huán)形鼓泡方案時(shí)，則需要進(jìn)一步提高噴氣壓力。

圖11 給出了t=0.3,0.9,2.4,12.5 s 時(shí)，中間的排氣孔中心截面的氣相云圖。相應(yīng)的溫度分布與中心孔鼓泡類似，因此認(rèn)為不存在結(jié)冰的可能性。由圖可以看出排出的氦氣逐漸往上移動(dòng)，隨時(shí)間增加，氦氣在液相中擴(kuò)散，所占體積越來越大，對液相存在強(qiáng)烈的攪拌作用，這現(xiàn)象與實(shí)驗(yàn)一致。計(jì)算結(jié)果表明在t=7.5 s 左右氦氣到達(dá)氣液界面進(jìn)入氣相空間。

圖10 環(huán)形和中心孔鼓泡氣枕壓力隨時(shí)間的變化Fig.10 History of ullage pressure of one hole and ringlike holes schemes

為進(jìn)一步證實(shí)氦氣的擴(kuò)散過程，圖12 給出了t=0.3,0.6,0.9 s 時(shí)環(huán)形小孔上方水平平面上氣相體積含量云圖。由圖可清楚地看到氦氣在水平方向上的擴(kuò)散，大概在0.6 s 時(shí)，相鄰的云團(tuán)開始聚合。一個(gè)問題是，為什么中心鼓泡時(shí)氦氣沒有擴(kuò)散，環(huán)形鼓泡時(shí)出現(xiàn)了嚴(yán)重的擴(kuò)散現(xiàn)象。定性的解釋是中心孔鼓泡時(shí)雖然噴出孔面積較大，但氣液之間的接觸面積要小得多，相同體積流量時(shí)噴出速度小，對貯箱煤油的流動(dòng)攪拌局限在小孔上方的圓柱形體積內(nèi)。但是環(huán)形鼓泡時(shí)，小孔總面積要小得多，噴出速度相對更快，氣液接觸面積也大得多，氣體更容易被液體攜帶，且相鄰小孔對流場的影響區(qū)域存在重疊，加劇了氦氣的擴(kuò)散。

4 結(jié) 論

本文基于Euler-Euler 兩相數(shù)值模型和SST k-ω湍流模型，對冷氦氣在煤油貯箱中的兩種鼓泡結(jié)構(gòu)中心孔式和環(huán)形小孔式進(jìn)行了數(shù)值模擬。模型與氦氣在水中的環(huán)形小孔鼓泡實(shí)驗(yàn)結(jié)果進(jìn)行了定性驗(yàn)證。對比了煤油和水鼓泡時(shí)氣枕壓力的變化。基于數(shù)值結(jié)果對煤油結(jié)冰可能性進(jìn)行了理論分析。得出主要結(jié)論如下。

（1）構(gòu)建的數(shù)值模型定性準(zhǔn)確地捕捉到了兩種鼓泡方式的流動(dòng)和傳熱特性。

（2）數(shù)值結(jié)果表明低溫氦氣進(jìn)入液體后在有限行程內(nèi)馬上升溫到接近煤油溫度，理論分析表明固化可能性微乎其微，實(shí)驗(yàn)觀察證實(shí)了這個(gè)結(jié)論。

（3）中心孔鼓泡中氦氣幾乎沒有擴(kuò)散，在噴射孔上方至液面下形成半徑比中心孔稍大的氣柱。但是環(huán)形小孔鼓泡方案數(shù)值和實(shí)驗(yàn)都觀察到了嚴(yán)重的氦氣擴(kuò)散現(xiàn)象，不滿足煤油排出要求。

圖11 t=0.3,0.9,2.4,12.5 s時(shí)中間的排氣孔中心截面的氣相云圖Fig.11 Phase contour of middle section of tank at t=0.3,0.9,2.4,12.5 s

圖12 t=0.3,0.6,0.9 s時(shí)環(huán)形小孔上方水平平面上氣相體積含量云圖Fig.12 Gas fraction contour of horizontal plane above ringlike holes at t=0.3,0.6,0.9 s

（4）由于氦氣和煤油換熱更充分，相同條件下煤油儲罐氣枕壓力比水儲罐氣枕壓力上升更快，但差別不大，可以用水代替煤油進(jìn)行冷氦氣鼓泡實(shí)驗(yàn)。