臥式液氫貯罐內(nèi)溫度分層數(shù)值模擬

王天祥，陳 虹，雷 剛，李愛華

(總裝備部工程設(shè)計(jì)研究總院，北京100028)

1 引言

隨著我國航天技術(shù)的發(fā)展，地面低溫推進(jìn)劑的貯存量越來越大，低溫真空絕熱貯罐容積也相應(yīng)增大。在低溫推進(jìn)劑大容量貯存過程中，由于貯罐的絕熱性能較高，進(jìn)入貯罐內(nèi)部的熱量較少，不能完全加熱貯罐內(nèi)的低溫推進(jìn)劑，使得在內(nèi)部近壁面區(qū)域的流體溫度升高，在浮升力的作用下沿壁面向上運(yùn)動，而內(nèi)部核心流體的溫度基本不變，從而在貯罐高度方向形成溫度梯度，即溫度分層現(xiàn)象.研究表明，溫度分層現(xiàn)象不但會導(dǎo)致貯罐內(nèi)推進(jìn)劑的上面層形成一定厚度的熱液層，影響推進(jìn)劑的輸出品質(zhì)，而且會惡化容器(貯罐)的使用性能，形成急劇的氣枕壓力上升［1］。

國內(nèi)外對低溫容器內(nèi)溫度分層的研究進(jìn)年來也有相關(guān)報(bào)道。文獻(xiàn)［2］采用渦流函數(shù)法對低溫流體液氫熱分層進(jìn)行了數(shù)值模擬，同時擬合出表面溫度對熱流密度、液體填充率及時間的函數(shù)表達(dá)式。文獻(xiàn)［3］利用邊界層積分方法對火箭液氫推進(jìn)劑的熱分層進(jìn)行了分析，并選用水作為研究工質(zhì)進(jìn)行了模擬實(shí)驗(yàn)。文獻(xiàn)［4］采用試驗(yàn)方法對大型容器內(nèi)液氫熱分層現(xiàn)象進(jìn)行了研究，并分析了振蕩與晃動對液氫熱分層現(xiàn)象的影響。文獻(xiàn)［5］通過對方型透明容器采用壁面加熱的可視化試驗(yàn)方法，進(jìn)行了水的熱分層現(xiàn)象實(shí)驗(yàn)研究，并且采用數(shù)值模擬的方法對圓柱型筒體內(nèi)液氫的熱分層規(guī)律進(jìn)行了分析，兩種結(jié)果對比得到了相似的結(jié)論。文獻(xiàn)［6］等針對兩個典型的液化天然氣渦旋事故，分析了分層和渦旋現(xiàn)象的機(jī)理，并且對雙向?qū)α鲾U(kuò)散模型、Bates-Morrison模型和四階段模型等進(jìn)行了介紹。文獻(xiàn)［7］分析了液化天然氣貯槽內(nèi)分層現(xiàn)象形成的原因、分層穩(wěn)定性和破壞性的機(jī)理以及分層的存在對各個液體層之間的傳熱和傳質(zhì)的影響，并給出了此種情況下的熱流率和質(zhì)流率的計(jì)算關(guān)聯(lián)式。文獻(xiàn)［8］采用大渦模擬方法對港口取水系統(tǒng)進(jìn)行了數(shù)值計(jì)算，得到了水溫與流速的分布。文獻(xiàn)［9，10］用FLUENT對火箭液氧豎直貯箱內(nèi)液氧溫度場和速度場進(jìn)行了數(shù)值模擬，得到了火箭液氧貯箱內(nèi)溫度與速度分布。文獻(xiàn)［11］采用CFD技術(shù)，對不同氣枕壓力下液氫貯箱內(nèi)部物理場進(jìn)行數(shù)值模擬。研究表明:氣枕壓力會對貯箱內(nèi)部溫度場與速度場產(chǎn)生重要影響。根據(jù)現(xiàn)有文獻(xiàn)資料，目前大多數(shù)研究主要集中在立式圓柱形容器的熱分層方面，對大容量臥式貯罐內(nèi)低溫推進(jìn)劑(尤其是液氫)溫度分層的研究很少。同時，低溫推進(jìn)劑臥式貯罐結(jié)構(gòu)龐大，加之液氫/液氧自身的物性特殊，進(jìn)行可視化等實(shí)驗(yàn)研究非常困難。

因此，本文采用CFD技術(shù)通過對大型臥式液氫貯罐內(nèi)溫度場的數(shù)值模擬，揭示貯罐內(nèi)部壓力場、溫度場的分布規(guī)律，分析了液氫熱液層的形成過程及原因，為發(fā)射場大規(guī)模低溫推進(jìn)劑的安全貯存和確保輸出推進(jìn)劑品質(zhì)提供了重要的理論支撐。

2 模型建立與邊界條件

液氫貯罐由水平圓柱體及左右兩橢圓型封頭組成，容積300 m3，采用內(nèi)外罐結(jié)構(gòu)，內(nèi)罐光滑，絕熱方式為高真空絕熱。由于時間跨度長，計(jì)算量巨大。由于儲氫罐是橫向放置，相對于垂直方向，水平方向的傳熱與流動可以忽略，故可以簡化為兩維問題。計(jì)算中采用二維結(jié)構(gòu)化網(wǎng)格，左右壁面用邊界層網(wǎng)格局部加密，計(jì)算網(wǎng)格如圖1所示。

圖1 模型網(wǎng)格

由于液氫貯罐內(nèi)存在自由相界面，計(jì)算中選用VOF模型和二維非穩(wěn)態(tài)粘性流體k－ε模型相結(jié)合的計(jì)算模型，采用熱邊界條件、標(biāo)準(zhǔn)壁面函數(shù)和基于壓力的計(jì)算方法。另外，由于液氫的特殊屬性，液氫密度采用Boussineqs模型，編寫了氫氣的密度、比熱、導(dǎo)熱系數(shù)等自定義函數(shù)。同時，針對FLUENT軟件自帶模型在解決相變問題時應(yīng)用范圍狹窄、模型也不夠完善的缺點(diǎn)，為了對液氫貯罐內(nèi)的物理場進(jìn)行較為準(zhǔn)確的模擬，編寫了相變模型的自定義程序。

描述貯罐內(nèi)部流動及傳熱的控制方程為:

式中 φ為廣義標(biāo)量，可以表示速度u，v，湍流脈動動能k，湍流脈動動能耗散率ε，溫度T。αp為氣態(tài)氫，αq為液態(tài)氫。ρ為密度，Γφ、Sφ分別表示擴(kuò)散項(xiàng)和源項(xiàng)。

由于外部熱侵使得貯罐內(nèi)部會發(fā)生氣液兩相的質(zhì)量和能量的變化，因此，為了較為真實(shí)的模擬貯罐內(nèi)部的物理現(xiàn)象，本模型通過定義連續(xù)性方程(1)中的和和能量方程中的源相來模擬液相蒸發(fā)過程，采用赫茲一努森方程，將其編寫成程序后外掛于Fluent進(jìn)行計(jì)算。

本研究的初始條件和邊界條件為:初始?xì)庹韷毫? atm，氣枕溫度為20.2 K，液體溫度為20 K，液氫填充率0.5;邊界為壁面邊界條件，輸入熱量為1.1 w/m2。

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

由于貯罐受到外部環(huán)境的熱侵，導(dǎo)致貯罐內(nèi)的液氫汽化，從而使得氣枕壓力升高。對于貯罐內(nèi)壓力的變化可以通過簡化理論給出計(jì)算值。假定氫氣和液氫在各自部分充分混合，溫度均勻，則可通過傳入貯罐內(nèi)的熱量，推導(dǎo)出能夠蒸發(fā)出的液氫量，從而根據(jù)氫氣狀態(tài)方程，計(jì)算出貯罐內(nèi)壓力變化，如圖2所示。圖2也同時給出了CFD計(jì)算的壓力隨時間變化率。從圖中可以看出，隨著貯存時間的增加，貯罐內(nèi)氣枕壓力逐漸上升，理論計(jì)算得到的壓力增長速率要小于CFD計(jì)算結(jié)果，這是由于在理論計(jì)算時沒有考慮貯罐內(nèi)液體的溫度分層，壁面漏熱被認(rèn)為是均勻的加熱整個液體，所以蒸發(fā)的液氫量相對于考慮溫度分層時要小。因此，在靜置貯存時，由于溫度分層會導(dǎo)致液氫貯罐內(nèi)的壓力急劇上升。同時，CFD計(jì)算結(jié)果顯示，壓力先開始有下降過程，是因?yàn)橛?jì)算時的初始值液氫未達(dá)到飽和，故氫氣向液氫傳熱導(dǎo)致氫氣溫度下降從而有壓力降低過程。

圖2 數(shù)值模擬與理論計(jì)算氣枕壓力比較

圖3 為靜置臥式貯罐橫截面上的流體流動速度矢量圖，從圖中可以看出，貯罐內(nèi)部流動是一個不穩(wěn)定狀態(tài)，渦旋的位置也在不斷變化。隨著貯存時間的增加，液相近壁處的邊界層內(nèi)熱液體在浮升力的作用下沿壁面向上運(yùn)動，積聚于氣液界面，而氣液界面下的冷流體向下運(yùn)動，形成了明顯的自然對流。由于受貯罐圓形橫截面的影響，熱液層內(nèi)液體質(zhì)點(diǎn)沿弧形軌跡運(yùn)動時，每一個點(diǎn)的受力情況與豎直貯罐內(nèi)的情況截然不同，從而導(dǎo)致冷液體在向下運(yùn)動時不完全是沿重力方向，其流線在不同時間具有不同的形狀，部分地方出現(xiàn)了更為復(fù)雜的渦旋流動。臥式貯罐內(nèi)的這種復(fù)雜流動，在一定程度上對液體的溫度分層起到了抑制作用，減小了容器壓力急劇上升的危險性，但對貯罐內(nèi)液體的輸出品質(zhì)(溫度)有一定的不利影響，因此，在液氫加注或輸送之前需要首先對貯罐泄壓，以保證液氫輸出品質(zhì)。

圖4為靜置臥式貯罐內(nèi)流體溫度分布圖，從圖中可以看出，臥式貯罐的邊界層覆蓋了整個與液體接觸的貯罐內(nèi)壁，由底部向液體表面厚度逐漸增厚。隨著貯存時間的增加，貯罐內(nèi)豎直方向上液體形成了較為明顯的溫度分層。在貯存過程中，通過貯罐內(nèi)壁面附近熱邊界向氣液界面不斷輸送熱液體，使得氣液界面下的熱液層不斷增厚，同時，由于貯罐內(nèi)液體的不穩(wěn)定流動，使得貯罐內(nèi)的主體液體的溫度也有少許增加。由于氣液界面一直處于動態(tài)平衡過程，其溫度也基本處于氣枕壓力下的飽和溫度，這可以通過圖2和圖4看出。

圖5是靜置臥式貯罐橫截面豎直方向上溫度隨時間變化曲線圖。隨著貯存時間的增加，氣液表面的熱液層厚度(氣液界面溫度與主體液體溫度之間的區(qū)域)逐漸增加，熱液層內(nèi)存在著一定的溫度梯度，其分布接近于高斯概率曲線，且氣相區(qū)也存在著非線性的溫度分布。在靜置貯存15小時后，熱液層的最高溫度與主體液體溫度相差1 K左右，而氣相最高溫度出現(xiàn)在頂部，約60 K左右。

圖3 靜置臥式貯罐內(nèi)流體速度矢量圖

圖4 靜置臥式貯罐內(nèi)流體溫度分布圖

圖5 靜置臥式貯罐橫截面豎直方向上溫度分布曲線

4 結(jié)論

(1)靜置臥式貯罐罐壁與液氫的換熱為帶邊界層的自然對流換熱，邊界層覆蓋整個與液氫接觸的貯罐內(nèi)壁面。

(2)靜置臥式貯罐內(nèi)的液氫沿豎直方向上出現(xiàn)溫度分層，其溫度分布近似于高斯概率曲線，最高溫度為氣液界面處的氣枕壓力下的飽和溫度，最低溫度為主體液氫的溫度。氣相區(qū)域存在著非線性的溫度分布，最高溫度出現(xiàn)在貯罐頂部。

(3)靜置貯罐內(nèi)的液氫處于不穩(wěn)定流動狀態(tài)，在局部地方有渦旋流動。臥式貯罐內(nèi)的這種復(fù)雜流動，使得貯罐內(nèi)的主體液體的溫度也有少許增加，在一定程度上對液體的溫度分層起到了抑制作用，減小了容器壓力急劇上升的危險性，但對貯罐內(nèi)液體的輸出品質(zhì)(溫度)有一定的不利影響，因此，在液氫加注或輸送之前需要首先對貯罐泄壓，以保證液氫輸出品質(zhì)。

［1］尼瓦費(fèi)林，亞波布拉諾夫.液體低溫系統(tǒng)［M］.北京:低溫工程編輯部，1993.

［2］Tanyun Z，Zhongp in H，Li S.Numerical simulation of thermal stratification in liquid hydrogen［J］.Advances in Cryogenic Engineering，1996，41:155 ～161.

［3］Tatom JW，BrowW H，KnightL H，et al.Analysis of thermal stratification of LH2 in rocket propellant tank［J］.Advances in Cryogenic Engineering，1964，9:265 ～272.

［4］Bailey T E，F(xiàn)earm R F.Analytical and experimental determination of LH2 temperature stratification［J］.Advances in Cryogenic Engineering，1964，9:254～264.

［5］DAS S P，CHAKRABORTY，DUTTA P.Studies on thermal stratification phenomenon in LH2 storage vessel［J］.Heat transfer engineering，2004，25(4):54～66.

［6］林文勝，顧安忠，李品友.液化天然氣的分層與漩渦研究進(jìn)展［J］.真空與低溫，2000，6(3):125～132.

［7］程 棟，顧安忠.液化天然氣的貯存分層現(xiàn)象［J］.深冷技術(shù)，1997(1):13～15.

［8］江春波，馬強(qiáng)，付清潭.二維溫度分層流的數(shù)值模擬［J］.水力發(fā)電，2003(2):24～26.

［9］程向華，厲彥忠，陳二峰.火箭液氧貯箱熱分層現(xiàn)象數(shù)值模擬［J］.低溫工程，2008(2):10～13.

［10］程向華，厲彥忠，陳二峰，等.新型運(yùn)載火箭射前預(yù)冷液氧貯箱熱分層的數(shù)值研究［J］.西安交通大學(xué)學(xué)報(bào)，2008(9):1132～1136.

［11］王磊，厲彥忠，程向華.氣枕壓力對液氫貯箱熱分層的影響規(guī)律［J］.低溫工程，2009(6):18～22.