貯箱內低溫推進劑汽化過程的CFD數值仿真

陳 亮 梁國柱 鄧新宇 胡 煒

(北京航空航天大學 宇航學院，北京100191)(北京宇航系統工程研究所，北京100076)

液體運載火箭貯箱在加注低溫推進劑后至發射前一直處于地面停放狀態，貯箱壁面漏熱會使推進劑的溫度升高并發生汽化，減少推進劑的可用量，并帶來安全隱患，因此需要對貯箱內傳熱傳質過程，尤其是汽化過程進行研究.國內外對貯箱內低溫推進劑的傳熱傳質過程的仿真采用的方法主要有集總參數法和CFD(Computational Fluid Dynamic)方法.文獻［1］對增壓貯箱進行了集總參數建模，將貯箱內部分為增壓氣體、液態推進劑主體部分和氣液界面3部分，根據熱力學平衡原理建立了相變模型.文獻［2］采用集總參數法針對液氧貯箱的增壓過程建立了貯箱內氣相空間的集總參數模型，認為氣枕與液氧表面間的傳質過程是由液氧的飽和蒸氣壓與氧氣分壓之差決定的.集總參數法在對貯箱的研究中應用比較廣泛，但不能反應貯箱中各項參數的空間分布情況.本文研究的貯箱尺寸較大且溫度分布對相變有較大影響，因而不適用于集總參數法.

國內外采用CFD方法的研究相對較少.文獻

［3］使用有限體積法研究通入熱增壓氣體后液氧貯箱內的物理變化，但忽略了相變過程，文獻［4］在此基礎上加入了相變模型，認為引起液氧相變的主要原因是增壓氣體對液氧的加熱作用，氣體向推進劑的傳熱是由自然對流引起的.參考文獻［5－6］對貯箱內液體主體區域采用二維的CFD方法建模并與氣枕的集總參數法模型耦合，計算得到的增壓速率與采用集總方法的熱力學分析結果相符.文獻［7］采用CFD方法建立了液氧貯箱的2維模型，根據熱力學平衡原理建立了相變模型.這些研究多采用類似文獻［4］的相變模型，氣液兩相之間換熱系數需根據流體的運動狀態求解，但在實際情況中，液態推進劑中可能出現大量“氣泡”，流動狀態復雜，不易求解.本文基于實際問題，根據基本的熱力學平衡原理構造相變模型，采用CFD方法對貯箱中的物理過程進行仿真.

1 物理過程分析

本文研究的推進劑貯箱為如圖1的高約7m、直徑約3 m的液氫貯箱，氫箱與氧箱共底，氫箱后底和氧箱前底間為真空隔熱層，貯箱內壁采用鋁合金，鋁合金壁面外包覆絕熱層.貯箱處于地面停放狀態，柱段外壁面直接與外界環境接觸，兩者之間的換熱過程為強制對流過程，受外界溫度和風速影響.氫箱前底外為火箭艙室，艙內充填氮氣，前底與艙室內氣體的換熱過程可視為自然對流過程.氫箱頂部開口，通過管路與外界大氣相通.外界環境以及氧箱(液氧溫度明顯高于液氫)通過壁面漏熱加熱貯箱中處于飽和狀態的液氫，使液氫發生汽化，這一過程屬于非均勻大容器沸騰［8］.在貯箱內部發生的物理過程主要包括推進劑的自然對流過程、推進劑與貯箱內壁面的對流換熱過程、貯箱壁的導熱過程以及推進劑的相變過程.其中，相變過程包括液體推進劑主體和氣枕間的界面上的蒸發和冷凝現象以及液體推進劑主體內部(主要是壁面附近)的汽化現象.

圖1 貯箱示意圖

2 數學模型

為滿足仿真需要，貯箱的數學模型除了包括多相流問題的通用控制方程，還需要建立推進劑相變模型、真實氣體模型和固體壁面導熱模型.

2.1 流體控制方程

控制方程中考慮流體粘性影響，忽略相變引起的動量變化.

連續方程:

式中，ρ為流體密度;v為速度矢量;Sm為質量源項;表示氣液兩相間傳質過程引起的單位體積的質量交換率.

動量方程:

能量方程:

式中，E為能量;keff為有效熱傳導系數為有效應力張量;Sh為能量源項，表示兩相間傳質過程引起的單位體積的能量交換率.

多相流方程

多相流模型選用可以追蹤相界面的VOF模型［5］，容積比率的連續方程為

式中，αq，ρq和vq分別為第q相的體積比率、密度和速度矢量;Sαq為第q相的質量源項.

2.2 推進劑相變模型

低溫推進劑的相變過程主要出現在加熱面(即貯箱內壁面)和氣液界面上(包括液體推進劑內部的氣泡界面和液體推進劑主體與氣枕間的界面).由于貯箱內推進劑近似處于飽和狀態，且貯箱壓力保持穩定，因此可假設液態推進劑及推進劑蒸汽處于熱力學平衡，則氣液界面上應滿足力平衡條件pg=pf，即氣液兩相壓力相等(不考慮張力影響);熱平衡條件Tg=Tf，即氣液兩相溫度相等;相平衡條件pg=pf=psat(T)，即氣液界面上壓力等于對應溫度下的飽和壓力［8］.在氣液界面上，當推進劑(i表示液氫或者液氧等低溫推進劑)的蒸汽分壓pig大于氣相區域溫度對應的飽和蒸汽壓psat，i時，發生冷凝，此時氣液界面上發生的質量交換和熱交換［1］分別為

式中，Δt為時間步長;Vig為氣相體積;Mi為推進劑摩爾質量;Til為液態推進劑溫度;R為氣體常數;Δhi為汽化潛熱.

當pig小于psat，i時，發生汽化.此時氣液界面上發生的質量交換率和能量交換率分別為

在貯箱壁面附近，壁面漏熱會引發推進劑的汽化.在發生汽化前，液態推進劑主體內部在內壁面附近沒有氣體，因此前述相變模型不適用.假設壁面附近處于熱力學平衡狀態，根據相平衡原理應滿足Tf=Tsat(p)，即液體溫度等于對應壓力下的飽和溫度，根據能量守恒的原理可得壁面附近的汽化過程的質量交換率和能量交換率分別為

式中，ρl為液態推進劑密度;Tsat，i為飽和溫度;cpi為液態推進劑的比熱容.

連續方程的質量源項Sm由式(5)、式(7)、式(9)決定;能量方程中的能量源項Sh由式(6)、式(8)、式(10)決定.由于相變模型建立在熱力學平衡假設的基礎之上，因此可用于近似處于飽和狀態且壓力溫度都不發生顯著變化的系統中，而對于過熱沸騰問題或壓力溫度發生劇烈變化的系統不適用.

2.3 其他模型

由于氫氣在低溫時的性質與理想氣體有一定差異，故根據試驗數據建立真實氣體模型如下:

式中，Z為壓縮因子，通過對試驗數據擬合得到，Z隨溫度升高而增大，當壓力為0.08 MPa，溫度為20 K 時，Z 約為0.9.

貯箱固體壁面的傳熱不僅包括垂直于壁面方向的傳熱過程，也包括沿壁面方向的傳熱，能量輸運方程如下:

式中，ρw為固體壁面的密度;hw為壁面材料的顯焓;kw為熱傳導系數;Tw為壁面的溫度.

3 計算方法

貯箱仿真模型通過Fluent建立.模型網格分為流體區域和固體壁面區域(包括鋁合金貯箱壁和絕熱層)，靠近固體壁面的流體區域網格進行適當加密.通過試算對比選定節點數為63 818的2維軸對稱網格，經驗證，不同劃分方式和節點數的網格的計算結果相近.

初始條件設置為推進劑處于飽和狀態.受重力影響，貯箱內不同深度的壓力不同，因此，將液氫的初始溫度設置為對應壓力下的飽和溫度.為簡化問題，貯箱固體壁面區域的初始溫度設置為與同一深度的液氫的溫度相同(實際情況中，在加注液氫前要對貯箱進行預冷，預冷后壁面溫度接近液氫溫度).

氫箱柱段和前底的外壁面設置為第3類傳熱邊界條件，壁面換熱系數根據對流運動狀態計算得到.環境壓力為0.08MPa，溫度為30℃，風速為3 m/s，計算得到柱段外壁面的強制對流換熱系數為9.008W/m2/K;氫箱前底外的火箭艙室內充填20℃的氮氣，自然對流換熱系數為 0.937～1.007 W/m2/K，隨外壁面溫度變化而變化.氫箱與氧箱的共底設置為第2類傳熱邊界條件，氧箱通過共底向氫箱的平均漏熱率根據試驗結果設置，為83 W/m2.湍流模型采用適合低雷諾數的RNG k-ε模型，壁面處理方式采用增強壁面函數.求解器選擇基于壓力的求解器，壓力插值方式選用body-force-weighted方案，速度壓力耦合方式選用壓力和速度耦合的coupled方案.相變模型通過UDF以質量源項和能量源項的形式加入到Fluent模型中.經多次試算比較，時間步長選定為0.04s，實際數值計算的連續性、動量、能量方程的殘差分別為 10－4，10－6，10－10量級.

4 計算結果與分析

貯箱的初始壓力、溫度和兩相分布如圖2所示.頂部壓力為0 Pa，底部為4957 Pa;頂部溫度為19.502 K，底部為19.697 K，溫度差為 0.195 K，液氫處于飽和狀態，呈明顯的溫度分層現象.隨著計算時間的推移，貯箱內推進劑出現對流運動和汽化現象.t=10.0 s和50.0 s時刻貯箱內溫度、局部速度、兩相和相變分布分別見圖3和圖4.由于側壁面和底面漏熱，靠近壁面的液氫溫度開始升高(如圖3a)，貯箱底部開始產生較明顯的汽化現象(如圖3c、圖3d，圖3d中黑色表示汽化，白色表示冷凝，灰色表示無相變).汽化產生的氣氫帶動液氫沿側壁面和底面向上流動(如圖3b).底部溫度較高的飽和液氫沿側壁面向上流動變為“過熱”，汽化加劇，上部液氫中出現大量氫氣(如圖4c);貯箱上部靠中間位置、溫度較低的液氫向下流動，變為“過冷”流體，汽化現象逐漸消失，并開始出現冷凝(如圖4d).

圖2 初始狀態的貯箱內各參數分布

圖3 t=10 s時刻貯箱內各參數情況分布

t=450 s時刻貯箱內溫度、局部速度、兩相和相變分布見圖5.

圖4 t=50 s時刻貯箱內各參數情況分布

圖5 t=450 s時刻貯箱內各參數情況分布

隨著壁面溫度的升高，外壁面與外界環境的換熱和內壁面與流體的換熱逐漸達到平衡，貯箱內液氫的溫度形成相對穩定的中間溫度較低，靠近側壁面和底面的區域溫度較高的“V型”溫度分布(如圖5a).貯箱內液氫的流動狀態也隨之趨于穩定，形成底部液氫沿側壁面向上流動，上方液氫從中間向下流動的循環流動(如圖5b).隨著溫度分布和流動狀態趨于穩定，貯箱內的汽化現象也趨于穩定，汽化現象主要集中在貯箱上部壁面附近(如圖5c，圖5d).

仿真中貯箱壁面初始溫度設置為與推進劑相同，壁面溫度受外界環境傳熱影響而上升，因此，外壁面換熱率減小，內壁面換熱率增大，換熱率變化曲線如圖6a所示，可見內外壁面的換熱率趨于相等，300s后基本穩定在約20000W.由于貯箱通大氣，箱內壓力接近環境大氣壓，因此貯箱出口流量近似等于液氫單位時間的汽化量(即汽化速率).貯箱出口流量變化曲線如圖6b所示，在450 s仿真中，液氫汽化速率在100 s以后開始趨于穩定(100 s以后平均值為23.4 g/s)，結合前文分析，液氫的相變與其溫度分布和流動狀態密切相關，汽化速率趨于穩定可以間接反映出隨著貯箱壁面傳熱過程逐漸穩定，液氫的溫度分布、流動狀態和相變情況也逐漸趨于穩定.

圖6 換熱量和出口流量變化曲線

假設貯箱壁面漏熱全部作用于相變過程，推進劑的汽化速率約為44 g/s，與仿真結果相比有出入，這一方面由于壁面漏熱中約有15%(與氣枕相接觸的內壁面對氣枕的換熱量約為3000 W)用于加熱氣體并被帶出貯箱;另一方面，由于貯箱尺寸較大而仿真時間相對較短，雖然推進劑汽化速率已經趨于穩定，但貯箱內溫度分布還繼續緩慢變化，一部分壁面漏熱加熱從上部流入底部變為過冷的液體而尚未直接影響到相變過程.因此，壁面漏熱雖是引起汽化的關鍵因素但并不能直接決定汽化速率，壁面漏熱對相變的影響和貯箱內各物理過程間的復雜的相互作用機理有待進一步研究.此外，相變模型建立在熱力學平衡假設基礎上，根據仿真，貯箱內推進劑的過冷度或過熱度在0.1 K以內，接近飽和狀態，但局部區域的不平衡狀態仍可能對仿真造成一定影響，因此相變模型也需要進一步改進完善，以解決非平衡狀態下相變過程的仿真問題.

5 結論

本文根據熱力學平衡原理建立了低溫推進劑的相變模型，通過對地面停放狀態下的液氫貯箱的仿真獲得了以下結論:

1)隨著貯箱壁面傳熱過程的穩定，貯箱內部熱力學過程趨于穩定.貯箱內液氫出現底部沿側壁面向上，上部由中間向下的循環流動;液氫溫度呈中間溫度低，側壁面和底面附近溫度高的“V型”分布.

2)受推進劑對流運動的影響，貯箱內出現汽化的區域主要集中在貯箱上部靠近內壁面處.在450 s的仿真中，液氫的汽化速率趨于穩定，100 s后平均值約為23.4 g/s.

3)影響液氫汽化速率的主要因素的壁面漏熱和推進劑對流運動.壁面漏熱使靠近壁面的液氫溫度升高發生汽化，還影響推進劑的對流過程，是影響汽化速率的關鍵因素.

References)

［1］Zilliac G，Arif Karabeyoglu M.Modeling of propellant tank pressurization［R］.AIAA-2005-3549，2005

［2］王贊社，顧兆林，馮詩愚，等.低溫推進劑貯箱增壓過程的傳熱傳質數學模擬［J］.低溫工程，2007，160(6):28－37 Wang Zanshe，Gu Zhaolin，Feng Shiyu，et al.Simulation of heat transfer and mass transfer in cryogenic propellant tank pressurization process［J］.Cryogenics，2007，160(6):28－37(in Chinese)

［3］Majumdar A，Steadman T.Numerical modeling of pressurization of a propellant tank［R］.AIAA-99-0879，1999

［4］代予東，趙紅軒.運用數學方法模擬推進劑貯箱增壓［J］.火箭推進，2003，29(3):34－40 Dai Yudong，Zhao Hongxuan.Simulation of propellant tank pressurization using mathematics method［J］.Journal of Rocket Propulsion，2003，29(3):34－40(in Chinese)

［5］Panzarella C H，Kassemi M.On the validity of purely thermodynamic descriptions of two-phase cryogenic fluid storage［J］.J Fluid Mech，2003，484:41－68

［6］Panzarella C H，Kassemi M.Self-pressurization of large spherical cryogenic tanks in space［J］.Journal of Spacecraft and Rockets，2005，42(2):299－308

［7］Cheng Xianghua，Li Yanzhong，Chen Erfeng，et al.Effect of return inlet on thermal stratication in a rocket tank［J］.Journal of Thermophysics and Heat Transfer，2010，24(1):112－122

［8］徐濟鋆.沸騰傳熱和氣液兩相流［M］.北京:原子能出版社，2001:215－220 Xu Jijun.Boiling heat transfer and gas-liquid two-phase flow［M］.Beijing:Atomic Energy Press，2001:215－220(in Chinese)