朱洪來， 趙立偉， 梁紅義， 魏健健

（1.北京控制工程研究所， 北京 100190； 2.浙江大學制冷與低溫研究所， 杭州 310027）

1 引言

為實現(xiàn)更高比沖、更大的運載能力，早期各國的運載火箭應用液氫／液氧低溫推進劑，并逐步利用低溫火箭末級或上面級開展短期空間飛行任務，成為了低溫空間飛行器的雛形。 液氫貯箱工作在空間復雜熱環(huán)境下，內(nèi)部壓力隨工作時間逐漸增加，因此壓力控制是液氫貯箱關鍵技術之一。液氫貯箱面對空間惡劣的熱環(huán)境，由于空間微重力造成內(nèi)部介質自然對流消失，加之液氫的熱導率低，從而液氫貯箱內(nèi)出現(xiàn)嚴重的熱分層現(xiàn)象［1］。 熱分層直接造成貯箱內(nèi)部氣液兩相局部過熱，引起液氫非必要氣化，從而引起低溫貯箱內(nèi)壓力異常變化。

1987 年，波音公司對貯箱內(nèi)液體由于混合引起的氣相液化以及液體消除分層問題進行了實驗研究［2］。 1993 年，NASA Levis 研究中心通過數(shù)值分析與實驗數(shù)據(jù)對比方法研究了液氫貯箱強制對流去分層問題［3］。 2016 年，NASA GRC 對液氫儲罐的蒸發(fā)率進行了2 次實驗［4］，實驗檢驗熱力學排氣與流體混合系統(tǒng)在地面液氫試驗中抑制氣相和液相熱分層的能力。 2017 年，Huang 等［5］、Wang 等［6］在室溫溫區(qū)TVS 模擬裝置采用R141b為氣液相變貯存介質，針對TVS 的混合模式和排氣模式提出了3 種運行控制策略，進行了貯箱內(nèi)液體溫度、熱分層以及排氣損失的研究。

為分析低溫貯箱內(nèi)的流場和溫度場的演化規(guī)律，多采用計算流體力學（Computational Fluid Dynamics，CFD）通用軟件Fluent 或CFX 等進行物理仿真分析開展了數(shù)值分析。 Grayson 等［7］基于Fluent 軟件開展了一系列液氫貯箱自增壓、循環(huán)去分層過程的CFD 數(shù)值研究。 考慮到箱體的軸對稱，數(shù)值模型為二維軸對稱幾何。 兩相流動多計算采用體積含量模型（Volume of Fluid，VOF），相比于對氣相和液相分別建立Navier-Stokes（NS）控制方程的Euler-Euler 模型，該模型為單流體模型，即假設氣液兩相流為一種混合流體，只對混合流體（Mixture fluid）建立相應的NS 守恒方程，其中相含量通過一個額外的相體積輸送方程計算得到。 由于是單流體方程，控制方程相對Euler-Euler 模型減少，因此計算速度更快，收斂性更好。缺點是該模型假設氣液之間沒有相對速度，即不考慮相間的相對速度。 但考慮到儲罐內(nèi)流動速度小，因此該假設引起的誤差可以忽略。

由于漏熱條件等因素影響貯箱的增壓速率緩慢，且受時間步長的限制，采用Fluent 軟件計算出在給定時間內(nèi)的溫度分布極為困難。 本文采用分區(qū)模型，將貯箱內(nèi)部空間分為過熱區(qū)、飽和區(qū)和過冷區(qū)三塊區(qū)域，離散化氣相／液相區(qū)質量和能量守恒方程和熱力學方程，確定初始條件，求得液氫貯箱內(nèi)熱分層與壓力分布。

2 理論分析

2.1 分區(qū)模型

液氫貯箱的三區(qū)模型示意圖如圖1 所示，貯箱為橢球封頭加圓柱段構型。 貯箱漏熱為第二類邊界條件的均勻熱流密度漏熱，其熱流密度沿貯箱軸對稱。 不考慮由于外界氣溫和絕熱材料不均勻性導致的漏熱隨時間的變化。

圖1 液氫貯箱的三區(qū)模型示意圖Fig.1 Three zone model of the liquid hydrogen tank

在貯箱內(nèi)部介質所占空間分為氣枕區(qū)、氣液界面區(qū)和液氫區(qū)3 塊區(qū)域，各區(qū)分別為一個單獨控制體。 基于貯箱和外界環(huán)境條件沿貯箱軸向的對稱性，氣枕區(qū)僅在貯箱軸向方向存在較大溫度梯度，但無明顯壓力梯度變化。 由于重力作用，不考慮氣枕與貯箱壁面有附著液體。 氣液界面區(qū)抽象為一層無質量的界面層，分別對氣枕區(qū)和液氫區(qū)進行傳熱和傳質。 貯箱內(nèi)流動速度小，暫不考慮該氣液之間相對速度。

2.2 控制方程

2.2.1 氣枕控制體

對于氣枕區(qū)域，通過質量守恒可得氣枕質量變化率為式（1）：

通過能量守恒，可得溫度變化率為式（3）：

式中，TU為氣枕溫度，其為貯箱高度和時間的函數(shù)，qenU為環(huán)境對氣枕區(qū)的熱流密度，qUI為氣枕區(qū)對界面的熱流密度，pU為氣枕區(qū)壓力，VU為氣枕體積，hgsat為飽和氣體焓，cvU為氣枕區(qū)定容比熱。

氣枕溫度為式（4）：

根據(jù)理想氣體狀態(tài)方程可知氣枕各參數(shù)關系為式（5）：

2.2.2 氣液界面控制體

氣液界面溫度TI可以通過修正Alabovskii 方程［8］得式（6）：

式中，K參數(shù)隨溫度變化，可以通過液氫參數(shù)表查取。

氣液界面氣化或冷凝質量可通過等能量階躍條件［9－10］求得，如式（7）所示。

式中，qIB為界面對液相區(qū)的熱流密度，hg和hl分別為氣氫和液氫的焓值。hg（TI，pU） 和hl（TI，pU） 分別為氣氫和液氫在界面溫度和氣枕壓力下的焓值。

2.2.3 液氫控制體

對于液氫區(qū)域，通過質量守恒可得式（8）：

液氫質量的初始條件為式（9）：

通過能量守恒可得式（10）：

式中，TB為液體溫度，qenB為環(huán)境對液體區(qū)的熱流密度，hl為液體焓，cpB為液氫區(qū)定壓比熱。 液體溫度的初始條件為式（11）：

液體體積變化率可以寫為式（12）：

式中，ρl為液氫密度。

液體體積的初始條件為式（13）：

2.2.4 補充條件

氣枕和液體充滿整個貯箱，因此可得式（14）：

式中，VT為貯箱容積。

氣枕區(qū)冷凝液質量流量可以表示為式（15）：

根據(jù)傳熱學，氣枕區(qū)對界面的熱流密度和界面對液相區(qū)的熱流密度可以寫為式（16）、（17）：

式中，AI為氣液交界面面積，kU和kB分別為氣相和液相的熱傳導系數(shù)。

2.3 溫度與壓力分布求解

模型通過REFPROP V9.1 軟件可獲得給定三區(qū)中任一區(qū)對應的溫度、壓力下焓、比熱、密度和傳熱系數(shù)等熱物理參數(shù)，用于方程迭代計算。利用MATLAB 軟件可編程求解液氫貯箱內(nèi)自增壓過程曲線、維持時間和內(nèi)部液位變化。 程序求解的流程圖如圖2 所示。

圖2 求解程序流程圖Fig.2 Solver flow chart

3 結果與討論

3.1 模型驗證

本文采用低比表面積的類球體液氫貯箱作為研究對象。 貯箱采用標準2 ∶1 橢球封頭，貯箱內(nèi)徑及貯箱總體高度均為2600 mm。 貯箱內(nèi)液氫初始溫度為20 K，外部漏熱為第二類邊界條件，熱流密度為0.83 W／m2。 圖3 為20%和80%兩種初始充滿率下，F(xiàn)luent 軟件和節(jié)點法的三區(qū)模型計算得到的結果對比。 結果可得，在Fluent 模擬的增壓過程的時間內(nèi)，三區(qū)模型與Fluent 的計算結果較為接近，對不同初始充滿率下的貯箱的增壓計算偏差不超過30 kPa。 因此，可以采用計算速率更快的三區(qū)模型來計算得出貯箱的增壓曲線。

圖3 2 種模型的增壓Fig.3 Pressure increase by fluent and three zone model

在模型中，氣液界面厚度dx應遠小于氣相區(qū)的高度（H－Hgs）。 目前Fluent 計算中，初始充滿率80%的算例在氣相溫度分層上形成近似線性的情況，因此這里采用初始充滿率80%在133 875 s 時刻的溫度分層數(shù)據(jù)來作為得出氣液界面厚度的參考值。 使用三區(qū)模型計算得到133 875 s 時刻的液相區(qū)平均溫度Tl，氣相區(qū)平均溫度Tg，氣相區(qū)飽和壓力Ps和溫度Ts，液位H，計算得出氣相區(qū)頂部溫度Tm，并與Fluent 計算結果進行比較，可得如圖4 所示的計算結果。

圖4 模型氣相區(qū)溫度分層Fig.4 Temperature stratification by fluent and three zone model

如圖4 所示，三區(qū)模型的計算結果和Fluent基本符合。 液相區(qū)溫度相比于Fluent 計算結果在某些高度部分存在0.2 K 左右的偏差。 位于罐底的溫度偏差是由于液相區(qū)均勻漏熱采用的第二類邊界條件造成的，位于氣液界面處的偏差則是局部傳質使得流場溫度不均勻導致。 關于氣相區(qū)的溫度分布，F(xiàn)luent 的計算結果顯示，計算時刻其還未發(fā)展為線性，因此和線性分布的假設存在偏差。

在圖4 的計算中， 氣液界面的厚度為27.8 mm，僅為貯箱高度0.5%，符合三區(qū)模型中兩相區(qū)為薄膜的假設。 與文獻［11］所提及的氫氣與液氫過熱度為30 K 的工程經(jīng)驗和液氫貯箱0.83 W／m2的漏熱條件，對80%初始充滿率的算例，氣液界面厚度與貯箱高度比也是基本一致的。

3.2 壓力分布

圖5 為使用一維節(jié)點法的三區(qū)模型計算得到的3 種不同初始充滿率貯箱內(nèi)部氣相區(qū)壓力以及貯箱內(nèi)液位隨時間的變化情況。 液位變化顯示，20%初始充滿率的貯箱變化趨勢與60%和80%初始充滿率貯箱的情況不同，隨時間逐漸減小。 說明低初始充滿率情況下液體膨脹速率要低于蒸發(fā)速率，而初始充滿率高的情況下則反之。 在壓力計算方面與圖5 的計算結果不同的是，在增壓前期（約3 h），80%初始充滿率的貯箱增壓速率最大，20%的增壓速率最??；然而在長期存儲（約200 h）中，80%初始充滿率的增壓速率在三者中為最小，20%初始充滿率的增壓速率最大。 對于20%、60%和80%初始充注率的工況，自增壓到0.5 MPa 所需的時間分別為125 h、240 h 和390 h。

圖5 增壓曲線和液位變化Fig.5 Pressure increase and level change curve

在初終壓以及漏熱都相同的情況下，增壓速率與所經(jīng)歷的時間呈負相關，經(jīng)歷的時間與不同初始充滿率下單位體積的總內(nèi)能變化量ΔU呈正相關。 因此，從仲氫T-v物性圖（圖6）中可得，3種不同初始充滿率的等比體積線代表貯箱從0.1 MPa 增壓至0.5 MPa 的熱力過程，而3 條線的初始和終態(tài)位置可查詢得到增壓過程的單位體積內(nèi)能變化量ΔU。 內(nèi)能變化量越大，意味著漏熱時間越長，即表示增壓速率越慢。 從圖中可得，80%初始充滿率的增壓速率最慢，20%初始充滿率的增壓速率最快。

圖6 仲氫的T－v 圖Fig.6 Para-Hydrogen T-v cure

3.3 溫度分布

基于三區(qū)模型可計算得出初始充滿率分別為20%、60%和80%情況下的溫度分布情況，如圖7所示。 由于氫的兩相物理參數(shù)的差異，圖中溫度分布曲線的轉折點出現(xiàn)在氣液界面上，轉折點前的液相溫度梯度明顯小于轉折點后的氣相溫度梯度。 相同漏熱情況下，填充率越高的氣相溫度梯度更大。 因此對于加注后未使用的液氫貯箱的溫度監(jiān)控更為重要。

圖7 溫度分層結果Fig.7 Temperature stratification results

4 結論

采用分區(qū)模型可以有效求解地面液氫貯箱內(nèi)氣液兩相的溫度和壓力分布情況，通過分區(qū)模型的求解過程可得如下結論：

1）不同初始充滿率的貯箱在給定漏熱條件下，三區(qū)模型對貯箱內(nèi)的壓力和液位進行計算結果與Fluent 的計算結果得到較好一致性；

2）在前人試驗和仿真結果基礎上建立的基于守恒定律和熱力學定理三區(qū)模型求解方法，在液氫貯箱內(nèi)溫度與壓力分布中保持計算精度條件下，大幅減低求解計算規(guī)模；

3） 三區(qū)模型可以計算不同填充率下、不同時刻液氫貯箱內(nèi)溫度與壓力分布，為貯箱的壓力控制和熱分層抑制提供了判定依據(jù)，同時為后續(xù)的熱力學排氣與流體混合過程提供了初始條件。