任金平， 于春柳， 任永平， 張 琪

(1.隴東學院化學化工學院， 慶陽 745000； 2.蘭州理工大學石油化工學院， 蘭州 730050)

1 引言

隨著探月工程、空間站建設和補給、火星探測計劃等項目的實施，人類對深空的探測不斷深入，低溫液體推進劑的貯存期也不再是幾天、幾周，而是幾個月甚至數(shù)年時間，不但要滿足運載火箭短時間發(fā)射階段的使用要求，而且要適應航天飛行器長時間在軌階段任務的需要[1]。 液氫、液氧、液態(tài)甲烷等低溫液體作為一種清潔的推進劑，含有人類空間生存必須的碳、氫、氧等物質，在未來深空探測中具有廣闊的應用前景。 由于低溫液體自身沸點溫度很低，在地面停放階段和在軌運行階段均會蒸發(fā)產(chǎn)生蒸發(fā)氣(Boil off Gas， BOG)，且隨著蒸發(fā)氣的不斷產(chǎn)生，容器內(nèi)壓力和溫度升高，對容器產(chǎn)生了不利影響[2－3]。 在太空微重力環(huán)境下，低溫液體的氣液界面長時間是相互摻混的，氣液存在狀態(tài)和位置也不確定[4－5]，為保持容器內(nèi)壓力穩(wěn)定，采用蒸發(fā)氣排出技術十分困難，而且蒸發(fā)氣排出一方面縮短了低溫液體貯存周期，另一方面造成燃料浪費和太空環(huán)境污染，同時對航天飛行器在軌運行姿態(tài)控制產(chǎn)生很大危害[2，6－7]。 因此，有必要對低溫液體的蒸發(fā)量進行控制，對容器內(nèi)蒸發(fā)氣進行有效處理。

從20 世紀60 年代開始，NASA 針對上述問題開展了深入的理論分析、數(shù)值模擬和試驗研究[7]。 國內(nèi)外針對上述問題的研究主要集中在先進絕熱結構及材料、熱力學排氣技術和主動制冷技術3 個方面。 主動制冷技術是基于低溫制冷機的各種形式的低溫液體蒸發(fā)氣再液化技術，包括低溫制冷機和冷凝器結合冷卻技術[8]、低溫制冷機和圓盤噴頭結合冷卻技術[9]、低溫制冷機和冷屏結合冷卻技術[10－12]、低溫制冷機和熱管結合冷卻技術[13－14]、低溫制冷機和導熱帶耦合冷卻技術[2]。 Ho 等[9，13－14]用CFD 軟件對采用低溫制冷機和熱管結合冷卻方式的液氫貯箱進行了瞬態(tài)熱分析，模擬研究了制冷機關閉和開啟2 種狀態(tài)液氫貯箱內(nèi)的溫度分布，還研究了泵噴管出口流體速度對貯箱內(nèi)流體熱分層的影響和泵噴管位置對貯箱內(nèi)液體混流和冷卻效果的影響；Kartuzova等[15]針對NASA 多功能氫試驗臺，數(shù)值模擬了低溫制冷機和噴桿結合冷卻方式的主動控壓過程。張磊等[16]采用Fluent 軟件對帶有浸沒噴射裝置的液氫零蒸發(fā)貯存系統(tǒng)進行模擬，預測貯箱內(nèi)溫度場分布，并分析了噴射裝置尺寸參數(shù)對溫度場的影響；王麗紅等[17]建立了微重力下低溫液體零蒸發(fā)貯存系統(tǒng)的3D 模型，對采用6 種不同換熱結構低溫貯箱的溫度場和速度場進行了模擬研究；賁勛等[18]以G-M 制冷機作為冷源，液氮為介質，通過換熱器對500 L 液氮貯存容器內(nèi)部輸入冷量，實驗結果表明對液氮貯存空間氣相區(qū)和液相區(qū)分別輸入冷量，均能抑制系統(tǒng)壓力上升趨勢，實現(xiàn)零蒸發(fā)貯存的目的，其中對于液相區(qū)輸入冷量效率更高，能夠在較短時間內(nèi)降低系統(tǒng)壓力。與上述其他形式低溫液體蒸發(fā)氣再液化技術相比較，基于低溫制冷機與導熱帶耦合作用的低溫液體蒸發(fā)氣再液化技術(即低溫制冷機和導熱帶結合冷卻技術)雖受制冷機工作時發(fā)熱的影響，但該技術工藝和涉及裝置簡單，控制參數(shù)少，在地面試驗容易實現(xiàn)，更便于在工業(yè)生產(chǎn)中推廣應用。

為深入研究低溫液體蒸發(fā)氣再液化過程中容器內(nèi)流體溫度場和速度場的變化情況，本文建立低溫制冷機與導熱帶耦合作用的低溫容器數(shù)值計算模型，采用數(shù)值模擬的方法，運用Ansys Workbench 軟件對地面停放階段制冷機開啟后低溫容器內(nèi)蒸發(fā)氣再液化過程進行模擬，分析容器內(nèi)流體溫度場和速度場的變化規(guī)律，為低溫液體蒸發(fā)氣再液化技術在實際工業(yè)生產(chǎn)中的應用提供理論依據(jù)。

2 模型描述

2.1 物理模型

圖1 為低溫制冷機和導熱帶耦合冷卻技術裝置示意圖。 制冷機直接安裝在低溫容器氣相空間頂部，導熱帶固定在制冷機冷頭上，與容器內(nèi)液相低溫液體和氣相蒸發(fā)氣接觸，對低溫液體冷卻和蒸發(fā)氣冷凝液化，蒸發(fā)氣再液化后的液滴在重力作用下自動流回容器。 該技術可選用的小型低溫制冷機有雙級斯特林制冷機、G-M 制冷機、脈管制冷機和SV 制冷機，一般制冷量需求大時選用斯特林制冷機，當要求振動極小且冷量又較小時則用G-M 制冷機或者SV 制冷機，脈管制冷機制冷量相對較小，要求容器容積和漏熱量都要小[19]。 制冷機運行可以通過容器內(nèi)壓力來控制，設定容器內(nèi)壓力范圍，當壓力達到設定上限值時制冷機開啟，蒸發(fā)氣開始液化，當壓力降低到下限值時制冷機關閉。 低溫容器采用雙層壁圓筒結構，由內(nèi)外容器和高真空變密度多層絕熱層組成，內(nèi)外容器材料S30408，絕熱層材料為隔熱紙＋鋁箔。 導熱帶作為制冷機冷量的傳輸介質，用熱阻小的固體導熱材料制成，被直接固定在制冷機冷頭上，作為熱耦合元件，要求導熱帶在低溫下具有導熱系數(shù)大、熱阻和密度小、抗振性能好、結構簡單、加工容易的特點。

圖1 低溫制冷機和導熱帶耦合冷卻技術裝置示意圖Fig.1 Schematic diagram of device for coupling effect of cryogenic refrigerator with thermal strap cooling technology

建立低溫制冷機與導熱帶耦合作用的低溫容器數(shù)值計算模型時，忽略其他與本文研究內(nèi)容無關的測量裝置及結構，通過轉化將低溫容器各部分漏熱量以熱流密度的形式加載到內(nèi)容器各個部位外表面，制冷機冷量同樣以熱流密度的形式加載到導熱帶上端表面，這樣簡化后的低溫容器數(shù)值模型由內(nèi)容器和導熱帶兩部分組成。 內(nèi)容器由筒體、上下封頭和徑管組成，容積為80 L，筒體為圓柱形，材料為S30408，內(nèi)徑為350 mm，厚度為2.5 mm，長度為665 mm；上下封頭均選用標準橢圓封頭，材料為S30408，內(nèi)徑為350 mm，總深度為113 mm，直邊段高度為25.5 mm，厚度為2.5 mm；頸管由不銹鋼管和聚氨酯泡沫塞組成，不銹鋼管子材料為0Cr18Ni9，外徑為100 mm，內(nèi)徑為95 mm，總長度為100 mm，聚氨酯泡沫塞直徑為95 mm，長度為100 mm；導熱帶采用紫銅箔逐層疊加焊接而成，長度為600 mm，寬度和厚度均為24 mm[20]。

2.2 數(shù)學模型

通過對地面停放階段制冷機開啟后容器內(nèi)部整個流體區(qū)域的求解，確定溫度場和速度場分布，需要構造的動力學控制方程除質量方程、動量方程、能量方程外，還需要構造與體積比率有關的多相流方程。 控制方程中考慮流體粘性影響，忽略相變引起的動量變化。 為了滿足數(shù)值模擬的需要，還需要建立系統(tǒng)內(nèi)低溫液體相變模型。

2.2.1 動力學控制方程

質量方程如式(1)所示。

式中，p為流體壓力(靜壓)，g為作用在微元體上的重力體積力，μl、μg分別為低溫液體和蒸發(fā)氣的動力粘度。

能量方程如式(3)所示。

式中，El、Eg分別為低溫液體和蒸發(fā)氣微團的總能量，kel、keg分別為低溫液體和蒸發(fā)氣的有效導熱系數(shù)，SE為能量源項，表示兩相間傳質過程引起的單位體積的能量交換率。

由于地面重力作用，容器內(nèi)存在氣液自由相界面，且相互間不存在互相穿插。 因此，選用軟件提供的可以追蹤相界面的流體體積函數(shù)多相流模型，與體積比率有關的質量方程為式(4):

式中，αi為i相的體積比率，ρi為i相的密度，ui為i相的速度矢量，Smi為i相的質量源項。

2.2.2 相變模型

容器內(nèi)流體相變主要出現(xiàn)在容器壁面、氣液界面和導熱帶周圍。 由于相界面附近的氣液兩相溫度、速度對氣液相間的傳熱和傳質影響明顯，考慮相間傳熱傳質機理。 由文獻[18]可知，對流體輸入冷量，在液相區(qū)產(chǎn)生的作用明顯快于氣相區(qū)，所以短時間內(nèi)蒸發(fā)氣量不會減少，假設容器內(nèi)壓力保持穩(wěn)定，低溫液體及蒸發(fā)氣始終處于熱力學平衡狀態(tài)，相界面處液體溫度和氣體溫度均等于氣枕壓力p對應的飽和溫度Ts。 計算過程中通過控制方程求解整個計算域的溫度分布，比較網(wǎng)格溫度T與Ts的相對大小判斷是否發(fā)生相變，若發(fā)生相變則通過在連續(xù)方程和能量方程中增加相應源項的方法來保證相間熱量與質量的轉移。 由Hertz Knudsen 公式與Clausius-Clapeyron 方程推導得相變過程數(shù)學描述關系式如下[21－22]。

若T≥Ts，液相蒸發(fā)，可得式(5):

式中，Sml、Smg均為相變時氣液界面質量傳遞速率，ζ與物性有關，為控制相變強弱的時間松弛因子，為了滿足收斂性需要，取值為0.1/s。

伴隨著質量的轉移，能量轉移速率如式(7)所示。

式中，γ為流體的氣化潛熱。

3 邊界條件和初始計算條件

參考前期研究低溫液體蒸發(fā)過程結果[20]:已知容器總漏熱量4.5 W，由徑管漏熱、上下封頭漏熱、筒體漏熱、徑向支撐漏熱、軸向支撐漏熱組成，徑管漏熱占容器總漏熱量的51%，徑向支撐漏熱很小忽略不計，由于組成徑管的泡沫塞絕熱效果非常好，忽略泡沫塞漏熱，但將徑管漏熱等效均勻設置在聚氨酯泡沫塞下表面，軸向支撐漏熱等效在下封頭表面。 通過轉化將各部分漏熱量以熱流密度的形式加載到模型各個部位表面，這樣泡沫塞下表面、上封頭表面、筒體表面、下封頭表面熱流密度分別為330.3 W/m2、1.5 W/m2、3.5 W/m2、1.5 W/m2。 根據(jù)容器總漏熱量及其他綜合因素考慮，低溫制冷機選用Sunpower 公司的CryoTel-GT 型斯特林制冷機，輸出功率10 W(恒定)，在計算中將導熱帶與制冷機接觸面簡化為恒定熱流表面，冷量以熱流密度(17 361.1 W/m2)的形式加載到該面上，導熱帶其余界面與液氮及蒸發(fā)氣邊界設置為耦合邊界。

兩相流體介質為液氮及其蒸發(fā)氣，假定液氮為不可壓縮流體，蒸發(fā)氣作Boussinesq 假設處理，容器初始充滿率為75%，液氮及其蒸發(fā)氣初始溫度設置為80 K，環(huán)境溫度為293 K。 由于容器內(nèi)流體溫度變化范圍不大，忽略導熱帶材料紫銅導熱率隨溫度的變化，對于紫銅箔導熱率、液氮、氮低溫蒸發(fā)氣均取80 K 對應的物性參數(shù)。 模擬過程中只考慮地面重力作用，不考慮表面張力和熱毛細效應對液氣界面的影響。

對再液化過程中的相變問題處理，采用C 語言編寫液化相變函數(shù)程序，通過用戶自定義函數(shù)(UDF)定義相變模型，蒸發(fā)氣與液氮兩相之間的傳質在體積分數(shù)方程中以源項的形式被指定，考慮潛熱的影響，將能量源項通過能量方程指定。

4 模擬結果及分析

4.1 溫度場隨時間變化過程

容器內(nèi)流體溫度場隨時間變化云圖如圖2 所示。 從圖中可以看出，內(nèi)容器內(nèi)流體溫度場沿導熱帶呈現(xiàn)對稱分布，在制冷機冷量的作用下，容器內(nèi)流體溫度開始降低，從容器頂部到底部溫度有明顯的分層現(xiàn)象，且液氮溫度降低速度快于蒸發(fā)氣，使得液氮主體區(qū)溫度普遍低于蒸發(fā)氣區(qū)，最低溫區(qū)出現(xiàn)在容器最底部，在20 s 時容器底部已開始出現(xiàn)過冷液氮，且隨著制冷機的持續(xù)作用，過冷液氮量逐漸增多。 分析可能產(chǎn)生該現(xiàn)象的原因:①液氮與蒸發(fā)氣相比較，傳熱系數(shù)明顯大于蒸發(fā)氣，使冷量在液氮主體區(qū)的傳輸比較快；②液氮與導熱帶的接觸效果好于蒸發(fā)氣，使液氮與導熱帶的接觸熱阻小于蒸發(fā)氣與導熱帶接觸熱阻；③吸收冷量后溫度降低的流體密度增大，在重力的作用下向容器底部運動，最后在容器底部形成最低溫區(qū)。 這表明對蒸發(fā)氣(氣相區(qū))和液氮(液相區(qū))分別輸入冷量，均能起到蒸發(fā)氣再液化的效果，其中對于液相區(qū)輸入冷量蒸發(fā)氣再液化過程效率更高，此結論與文獻[18]對該技術進行的實驗研究結果一致；同時還表明該再液化技術最先產(chǎn)生效果的是低溫容器底部液相區(qū)，然后逐步從下到上，最后作用于氣相區(qū)。

圖2 容器內(nèi)流體溫度場隨時間變化云圖Fig.2 Cloud diagram of temperature field changes of fluid in vessel at different times

從圖中也可以看出，在氣相蒸發(fā)氣區(qū)，5 s 時有明顯的溫度分層現(xiàn)象，20～120 s 蒸發(fā)氣區(qū)溫度趨于均勻，分層現(xiàn)象消失，基本恢復到初始溫度，說明在該時間段內(nèi)，蒸發(fā)氣對容器壁漏熱的吸收強于對導熱帶冷量的吸收，而且吸收的漏熱量大于吸收的冷量。 從5～120 s 的溫度云圖可以看出，紫銅箔導熱帶內(nèi)部冷量沿軸向傳輸速度比沿徑向快的多，因此在該時間段內(nèi)導熱帶周圍未形成明顯的溫度場分布。 從120～720 s 溫度云圖可以看出，導熱帶開始與周圍流體之間進行傳熱，在導熱帶周圍也開始逐漸形成沿徑向和軸向分布明顯的溫度場，而且沿這兩個方向分布的溫度場不斷地向周圍擴展，軸向的擴展速度比徑向的快，徑向液氮區(qū)擴展比蒸發(fā)氣區(qū)快，但在蒸發(fā)氣區(qū)和液氮主體區(qū)上部靠近容器壁處，受容器壁漏熱的影響，流體溫度與初始溫度相比較無明顯變化。 氣液界面之間的傳熱傳質對再液化過程流體溫度場的影響很小，幾乎未表現(xiàn)出來。

4.2 速度場隨時間變化過程

容器內(nèi)流體速度場隨時間變化云圖如圖3 所示。 從圖中可以看出，5～40 s 容器內(nèi)各處流體流動相對都比較劇烈，在導熱帶兩側形成了多個渦旋，并且這些渦旋沿導熱帶呈現(xiàn)對稱分布，分析原因，開始制冷機產(chǎn)生的冷量集中在容器中心軸線位置，而在容器壁處有環(huán)境漏熱，這使容器內(nèi)沿徑向從流體中心到容器壁之間形成了比較大的溫差，該溫差導致了上述結果。 流體速度大小分布呈現(xiàn)導熱帶下端兩側液氮流動速度最大，氣相區(qū)蒸發(fā)氣流動速度次之，容器底部液氮流動速度最小的規(guī)律；5 s時最大速度區(qū)域出現(xiàn)在氣液界面處和導熱帶靠近下端處，此時流場擾動最為劇烈，隨著時間的推移，容器內(nèi)流體速度按照從上到下的順序逐漸減小，同時流場擾動減弱。 在80 s 時只有導熱帶下端到下封頭內(nèi)表面處液氮有明顯變化的速度場分布，而導熱帶所處位置由于制冷機冷量的作用，流體吸收的熱量不足以維持渦旋流動，渦旋消失，流體溫度趨于均勻，流體流動減緩。 從120～720 s 速度場云圖可以看出，由于制冷機的持續(xù)作用，容器內(nèi)流體溫度一直處于均勻狀態(tài)，流體流動非常緩慢，只有下封頭內(nèi)表面處局部液氮有明顯的速度場分布。同樣，從圖中也可以看出，氣液界面之間的傳熱傳質對再液化過程流體速度場的影響很小。

圖3 容器內(nèi)流體速度場隨時間變化云圖Fig.3 Cloud diagram of velocity field changes of fluid in vessel at different times

5 結論

本文研究了基于低溫制冷機與導熱帶耦合作用的低溫液體蒸發(fā)氣再液化技術，建立了低溫容器計算模型，運用Ansys Workbench 軟件對制冷機開啟后低溫容器內(nèi)蒸發(fā)氣再液化過程進行了數(shù)值模擬，分析容器內(nèi)流體溫度場和速度場變化情況得到下述結論:

1)基于10 W 小型斯特林制冷機與紫銅箔導熱帶耦合作用的冷卻技術，能夠使低溫液體冷卻和蒸發(fā)氣冷凝再液化，實現(xiàn)對容器內(nèi)低溫液體蒸發(fā)量的控制。

2)對氣相區(qū)和液相區(qū)分別輸入冷量，均能起到制冷的效果，其中對于液相區(qū)輸入冷量蒸發(fā)氣再液化過程效率更高，在容器底部最先出現(xiàn)過冷液氮，且量逐漸增多。

3)開始流體流動比較劇烈，形成了沿導熱帶兩側對稱分布的多個渦旋，隨后渦旋逐漸消失，流動減緩，溫度趨于均勻。

4)制冷機的作用使容器內(nèi)沿徑向從流體中心到容器壁形成了比較大的溫差，導致容器內(nèi)流體開始流動比較劇烈，且形成了沿導熱帶兩側對稱分布的多個渦旋，隨著制冷機的持續(xù)作用，渦旋消失，流體流動減緩，溫度趨于均勻，只有下封頭內(nèi)表面處局部液氮有明顯變化的速度場分布。

5)紫銅箔導熱帶內(nèi)部冷量沿軸向傳輸?shù)乃俣缺妊貜较蚩欤瑲庖航缑鎸υ僖夯^程中流體溫度場和速度場分布的影響很小。