液氫貯箱停放過程中的力熱分析

羅天培，張 偉，李 茂，張家仙

(1.北京航天試驗技術研究所，北京 100074；2.北京市航天試驗技術與裝備工程技術研究中心，北京 100074)

0 引 言

運載火箭在完成低溫推進劑加注后，由于貯箱壁面的漏熱，推進劑的溫度會發生變化并發生氣化現象，這會導致推進劑的量不斷減少[1]。尤其遇到緊急情況需要推遲發射并長時間停放時，了解推進劑的溫度品質變化規律、不同液位時推進劑的蒸發率變化情況等，對于發射流程的制定至關重要。

早期對貯箱內推進劑和氣枕間的傳熱傳質分析一般采用集總參數法[2-3]，該方法具有較高的計算效率[4]，但需引入一系列的假設，無法全面模擬各物理量的空間分布情況。隨后，CFD技術以其高精度、流場信息全面、使用便捷等特性被廣泛地用于貯箱內的流場仿真。1990年，Hardy和Tomsik首次使用CFD技術分析了排液前不同參數對貯箱內溫度分布的影響，通過與試驗數據的對比，證明了CFD技術以及Flow-3D軟件在貯箱內流場仿真的有效性與準確性[5]。2015年，Wang和Li等人為了研究無排氣加注時貯箱內的熱力學變化情況，開發了一種可精細模擬貯箱固壁和流體間傳熱的CFD模型[6]，該模型可識別并模擬自然對流、核態沸騰、過渡沸騰以及膜態沸騰等不同傳熱模式，并適用于正常重力以及微重力。Fu和Sunden等人采用Lee相變模型模擬了微重力條件下貯箱的自增壓過程，文中詳細分析了重力、氣泡接觸角以及表面張力等作用的影響[7-8]。2016年，Liu和Li等人同樣利用Lee模型模擬了末級火箭貯箱受氣動加熱以及空間輻射的情況下的增壓性能以及熱分層現象，并開展了多層絕熱下貯箱的性能研究[9-10]。Kassemi和Kartuzova在模擬液氫貯箱自增壓的過程中，詳細分析了近平衡動力學模型中調整系數的作用以及氣、液間湍流環流的影響[11]。Wang和Ma等人利用CFD手段研究了貯箱金屬壁面設置內外兩層絕熱層對增壓氣體用量的影響[12]，結果表明當在貯箱內壁設置絕熱層時，可以顯著減少增壓氣體的用量。另外，除了單純地對貯箱內的流動及傳熱過程進行模擬以外，楊旦旦、李青等人還進行了貯箱流固耦合的計算，分析了液體晃動等情況對貯箱結構產生的影響[13-14]。

從國內外文獻來看，目前模擬氣液間相變的模型主要有兩類。第一類認為貯箱內的相變發生只發生在氣液交接面及固液交接面(貯箱壁與液體)上[15-16]，以熱力學平衡狀態作為相變判斷條件，當兩個交界面的熱力學平衡狀態被打破時即有相變發生；第二類以經典的Lee模型為代表[8,17]，以溫度為判斷條件，每迭代一步，對流場內全體網格進行判別，如果液體高于設定的相變溫度，則發生蒸發相變，相反，若氣體低于設定的相變溫度，則發生冷凝相變。顯然，第一類模型有其局限性，雖然從貯箱壁面的漏熱會先加熱固液交接面的液體，但隨著停放時間的加長，熱量會慢慢傳至流體內部區域，而該模型對于流體內部的氣化區域無法預測，隨著貯箱停放時間的延長，仿真精度勢必會越來越差。總體來說，利用該模型預測的蒸發量會偏小。第二類模型相對更為合理，其邏輯判斷由全場出發，但依然有其不足之處：其中相變溫度為定值，一般給為氣枕壓力下對應的飽和溫度，這在氣液交接面附近判斷較為合理，但由于重力作用，越接近容器底部的液體壓強越高，其對應的飽和溫度也越高，由氣液交接面的飽和溫度作為全流場相變的判斷依據，顯然會使預測的蒸發量偏大。

本文為探尋某型運載火箭液氫加注后貯箱內的力、熱變化規律，比較不同液位時貯箱的蒸發率情況，采用CFD分析方法對貯箱內部的流場進行求解，其中為了提高對相變過程的模擬精度，本文利用安托因方程對Lee模型進行了修正。以期尋找到長時間停放過程中最經濟液位(蒸發率最小)，為發射流程的制定以及緊急處置提供參考依據。

1 數學模型

為了模擬液氫的氣化過程，本文的計算需引入相變模型。引言中指出，Lee模型可以對流體內部的氣化區域進行模擬，更符合真實情況，因此本文選用了該模型。Lee相變模型[8]以貯箱內氣枕壓力對應的飽和溫度Tsat為判斷條件，每迭代一步，對流場內全體網格進行判別。當網格溫度Tcell≥Tsat時，液相蒸發，蒸發量表示如下：

(1)

當Tcell

(2)

伴隨著質量的轉移，能量轉移速率Sh表示如下：

Sh=-mvΔh

(3)

式中：mv、ml分別為氣相、液相質量轉移速率，kg/(m3·s)；Sh為能量轉移速率，J/(m3·s)；av為氣相體積分數；ρv、ρv分別為氣相和液相的密度，kg/m3；C為蒸發/液化系數，本文選取默認值0.1。

其中Tsat為定值，重力作用導致的容器底部壓強升高會使其對應的飽和溫度也升高，因此在容器高度較大或者推進劑密度較大時會使得蒸發量的預測偏差增大。為了考慮壓力和飽和溫度的對應關系，在此引入安托因(Antoine)方程進行修正，安托因方程表示如下：

lgP=A-B/(Tsat+D)

(4)

式中：P為壓強，Pa；A、B、D為物性常數，不同的推進劑對應于不同的值。

每次計算迭代中，先由全場網格的壓力值計算對應的飽和溫度，再分別對每個網格進行相應的相變判斷。

由于貯箱內始終存在著氣液相界面，因此兩相流模型選用適用于氣、液相間存在明顯分界面的VOF模型。VOF模型通過引進相體積分數實現對每一個計算單元相界面的追蹤，在每個計算單元內，所有相體積分數總和為1。引入VOF模型后，通過求解質量守恒、動量守恒、能量守恒三大方程獲得整個流場的信息。

此時，流體的質量守恒方程如下：

(5)

式中：aq、ρq、vq分別為第q相的體積分數、密度和速度矢量；mq為第q相的質量源項，對于氣相有mq=mv，對于液相有mq=ml。

動量方程為：

(6)

式中：μ動力黏性系數，Pa·s;g為重力加速度，m/s2。

對于VOF模型，兩相的能量方程是統一的，具體表示為：

(7)

式中：E為能量，J;keff為有效傳熱系數，W/(m2·K);teff為有效應力張量，Sh為能量源項，與式(3)保持一致。

為了封閉求解方程組，同時為了模擬貯箱內部的湍流效應，計算需引入湍流模型。按文獻[15]的方式分析結果，本文湍流模型采用適合于低雷諾數兩相流的可實現k-ε模型，其表達式如下：

Gk+Gb-ρε-YM+Sk

(8)

(9)

式中：k為湍動能，m2/s2；ε為湍流耗散率，m2/s3；μt為湍流黏性；Gk為由平均速度梯度生成的湍動能；Sk和Sε為源項。

2 仿真模型驗證

為了證明模型的準確性，利用NASA于1965年公開的一組液氫貯箱增壓排液過程記錄的試驗數據[18]開展了測試工作。

貯箱結構如圖1所示，箱體由圓柱筒體與上下橢球封頭組成，初始時刻箱內裝有溫度為25.6 K的液氫，增壓氣為氫氣，氣枕壓力1.1 Mpa，增壓排液過程中維持氣枕壓力不變。模型具體參數如表1所示。

圖1 貯箱模型Fig.1 Model of tank

計算過程中除了相變模型外，還編寫了自動增壓程序以維持氣枕壓力恒定，具體為：①當氣枕壓力超過設定的上限值時，降低增壓氣的入口流量；②當氣枕壓力低于設定的下限值時，增加增壓氣的入口流量；③當氣枕壓力在合理區間之內時，維持當前流量。

表1 貯箱主要參數Table 1 Main parameters of tank

文獻中分別給出了排液開始后90 s、178 s以及320 s的貯箱中線溫度分布，仿真計算后與試驗對比的情況如圖2所示，其中的曲線為計算結果，點為試驗獲得的數據。從圖中可見，3個時間點上仿真結果與試驗結果符合的均較好，3個時間點的平均仿真誤差分別為13.9%、9.9%及14.4%。由于本程序中相變發生的判據為溫度，而相變過程的吸熱及放熱也會直接影響溫度場，因此該算例的結果在一定程度上驗證了相變模型的準確性，后文用該模型開展的計算分析具備合理性。

圖2 不同時刻仿真與試驗的溫度場對比情況Fig.2 Comparison of temperature distribution between simulation and test at different time

3 貯箱停放過程模擬

3.1 計算方法及邊界條件

根據模型特點建立了二維軸對稱模型，壁面網格進行加密。由文獻[1]查知，氫箱柱段和上底的外壁面設置為第三類傳熱邊界條件，考慮環境溫度及風速，計算得到柱段外壁面的強制對流換熱系數為9.008 W/m2·K-1，氫箱上底外的火箭艙室內充填20 ℃的氮氣，自然對流換熱系數取為1 W/m2·K-1，氫箱與氧箱的共底設置為第二類傳熱邊界條件，氧箱通過共底向氫箱的平均漏熱率根據試驗結果設置，為83 W/m2。為了驗證網格無關性，分別選取總數分別為34022、62102、125837的三套網格開展計算。取貯箱中線的溫度分布作為對比，三套網格計算結果如圖3所示。從圖中可見，數量為34022網格的計算結果和其他兩組算例結果偏差較大，而62102及125837兩套網格的計算結果偏差很小，考慮到計算成本的節省，本文利用總數為62102的網格開展計算。

圖3 網格無關性測試Fig.3 Grid independence test

計算中氫氣密度采用理想氣體模型，為了考慮浮升力對液相內對流換熱的影響，液氫密度選擇Boussinesq模型，定壓比熱容Cp對溫度變化較為敏感，計算中采用線性擬和的方式給入。絕熱層導熱系數為0.02371 W/m·K，厚度給為20 mm，初始時刻氣枕壓力為0.1 Mpa，整個流場統一設為此壓強下對應的飽和溫度，隨著求解的進行，氣枕壓力和流場溫度會不斷趨近于實際情況，求解器選擇給予壓力的求解器，壓力插值方式選用body-force-weighted方案，體積分數項采用Geo-Reconstruct格式，速度壓力耦合方式選用壓力和速度耦合的coupled方案。經過多次試算比較，時間步長選定為0.02 s，實際數值計算的連續性、動量、能量方程的殘差分別為10-6、10-7、10-12量級。

氫箱總容積約為46.5 m3，本文共取了兩種典型的液位進行計算：①充填率50%，②充填至Ⅲ液位(44.8 m3)。依據文獻[1]所述，貯箱停放450 s后內外壁達到熱平衡，蒸發相變開始穩定，本文50%充填率算例計算時長達到了600 s，滿足了穩定相變過程發生對停放時間的要求，而對典型的Ⅲ液位連續計算時間達到24 h。

3.2 計算結果及分析

圖4給出50%充填率情況下貯箱內不同時刻的兩相流狀態。初始時刻(0 s)氣液交界面下為純液態的飽和液氫。在20 s時，從圖中可見貯箱壁面處液氫開始氣化，出現了成片的氣泡，并且隨著時間的推移，貯箱壁面的氣泡一直存在，這是由于貯箱壁面向內不斷漏熱的結果；在200 s時，隨著冷熱流體的不斷流動和熱交換，液體核心處也有液氫超過當地飽和溫度，開始氣化，同時氣液交界面下也開始出現氣泡；在450 s時，氣液兩相流狀態和200 s并無明顯的規律性變化，說明200 s時蒸發相變已進入穩定狀態。

圖5給出50%充填率40 s時刻貯箱內的矢量圖。從圖中可見，由于計算時考慮了浮升力的作用，壁面附近液體受熱后密度減小，裹挾著氣化后的氫氣向上移動，上升到氣液交界面時轉向貯箱中部移動，一周的熱流體運動至貯箱中線后對撞再向下移動，形成了一個大的旋渦，熱流體流動過程中會和貯箱內原有的冷流體不斷進行對流換熱，進而使貯箱內部液氫也開始出現氣化，與圖4中反映的物理過程相符。

圖6給出50%充填率40 s、200 s以及450 s時刻貯箱內的溫度分布云圖，初始時刻(0 s)貯箱內氣枕溫度設為當地大氣壓下對應的液氫飽和溫度，隨著時間的推移，貯箱氣枕空間的溫度場呈現了經典的分層分布現象。

圖7給出初始液氫加注到Ⅲ液位時，在距離貯箱頂部3.7 m的水平面上不同徑向位置溫度隨時間的變化，圖中0代表軸線位置，1.495接近貯箱壁面。從圖中可見，貯箱中推進劑不僅在水平高度上呈現溫度分層現象，在不同的徑向位置也呈現溫度分層的現象，其中越靠近貯箱壁面的推進劑溫度越高。從距離貯箱軸線1.495 m的位置，約前50 s溫度一直在升高，50 s開始溫度有個大幅下降，這是因為在50 s附近液氫發生了相變，由于氣化相變過程需要吸收大量的氣化潛熱，使得周圍液氫溫度反而下降，而50 s時其他三個位置相對于初始時刻(0 s)的液氫溫度有下降，這是由于貯箱內“旋渦”的作用，將冷流體卷動，通過對流換熱的作用降低了該處推進劑的溫度，隨后也開始出現類似1.495曲線的“上升—下降”過程，原因與1.495曲線類似，也是流體被加熱，隨后發生氣化降溫的過程，由于這三個位置在貯箱內部，相對于外壁附近的流體來說，該過程有一個約50 s的延遲。而在150 s后受熱、氣化過程區域穩定，各點溫度也開始變得平穩。

圖8給出貯箱從加注到Ⅲ液位開始停放24 h內蒸發率隨時間的變化，停放開始時蒸發率為1.2 m3/h，停放4 h左右，液位降至37 m3左右，此時蒸發率達到峰值，超過 2m3/h，在4 h至16 h，蒸發率慢慢減小，當停放至17 h附近時，蒸發率開始波動，此時貯箱內液氫剩余17 m3左右。

整個漏熱蒸發過程本質上為漏熱—推進劑溫度升高直到達到相變溫度—相變同時吸收氣化潛熱這樣一個過程，推進劑蒸發率的根本決定因素為氣液交界面傳熱量及固液交界面傳熱量。17～37 m3位于貯箱的柱段，此時氣液交界面面積一定，同時由于貯箱敞口，氣枕壓力基本不變，導致氣液交界面溫度基本不變，氣液交界面的換熱量也就基本相同。而停放過程中貯箱共底的換熱量也一致，這樣蒸發率的決定因素即為柱段的濕邊面積(固液交界面面積)，顯然柱段的幾何特征決定了濕邊面積與充填率成線形的關系，故物理過程完全類似，蒸發率也幾乎呈線性變化；而在17 m3前及37 m3后貯箱兩端的橢球形封頭為主要的幾何影響因素，在貯箱液位在17 m3左右時，貯箱的共底使得該區域的液體隨著充填率的增加，濕邊面積并未線性增加，但推進劑的整體熱容卻在增加，即從箱壁漏進來的熱量“更難”將推進劑加熱至其沸點之上，故整體的蒸發率呈現了一個短暫的平臺區和局部的極值，而在貯箱的充填率達到37 m3以上時，貯箱內的推進劑已注入到上封頭，此時氣液交界面面積在不斷減小，使得氣液交界面的傳熱量在減小，另外，除熱容的影響之外，封頭的漏熱系數也遠小于柱端(封頭1 W·m-2·K-1，柱段9.008 W·m-2·K-1)，使得隨著推進劑的增加，蒸發率反而在減小。

故若試圖使得貯箱在停放過程中保持較低的蒸發率，需將貯箱液位充填至37 m3以上或17 m3以下。

圖4 50%充填率貯箱內不同時刻的兩相流狀態Fig.4 The status of two phase flow in 50% filling rate at different time

圖7 同一水平面上不同徑向位置溫度隨時間的變化Fig.7 Temperature variation at different radial location in the same level

圖8 貯箱停放24 h過程中氣化率的變化Fig.8 Evaporation rate variation while parking in 24 h

4 結 論

本文利用CFD技術對某型運載火箭液氫貯箱的停放過程進行了模擬，其中相變模型采用了經安托因方程修正的Lee模型，得出的主要結論如下：

1)仿真模型對貯箱內的溫度場預測較為準確，3個時間點的平均誤差最大為14.4%，證明了采用安托因方程修正模型的合理性與有效性。

2)在貯箱停放一段時間(約450 s)后，在豎直方向與徑向均存在溫度分層的現象；液相內會形成大的漩渦，使得冷熱流體不斷進行熱交換，并導致貯箱內部的液氫也出現氣化。

3)貯箱停放期間蒸發率最大值超過2 m3/h，發生在停放4 h左右；而貯箱液位充填至37 m3以上或17 m3以下時蒸發率較低，最小值在12 m3/h左右，該數據可為發射流程的制定以及緊急處置提供參考依據。