飛行器用液氫貯箱低溫推進劑晃動熱力學響應分析

李 文,陳叔平,趙高逸,謝高峰

(1.蘭州理工大學石油化工學院,蘭州 730050; 2.北京航天發射技術研究所,北京 100076)

1 引言

隨著航空、航天技術的快速發展,飛行器的探測范圍和續航時間明顯增加,飛行器所攜帶液體推進劑的比重不斷增大,飛行器在姿態調整過程中貯箱內推進劑發生劇烈晃動,即充液容器在外部激勵作用下自由界面波動或容器強制運動所引起的流體整體運動現象[1-2]。而低溫推進劑密度大、粘度小,外部微小振動或晃動激勵即可引起罐內流體大幅非線性晃動,嚴重威脅飛行器系統穩定性和結構安全性,且該現象自20 世紀五六十年代就已引起國內外研究人員的廣泛關注[3-5]。

在一定漏熱環境下,外界擾動所產生的低溫推進劑晃動不僅會對貯箱產生附加力矩,還將引發貯箱內氣液相間熱力學不平衡等現象。劉展等[6-7]針對非等溫貯箱晃動動態響應進行了分析。Isaacson 等[8]針對圓柱形貯箱作簡諧運動和不規則運動時的液體晃動特性,定性分析了外部激勵與自由液面波形的關系。Ludwig 等[9]在圓柱形貯箱中用液氮進行了實驗,根據不同的晃動試驗,建立了Nusselt 數和Reynolds 數,以反映箱體壓降與晃動特性之間的關系。Das 等[10]采用FC-72 與HFE7000 為工質,實驗研究了在穩定區以及非穩定區,軸向激勵對箱體熱力特性的影響,研究指出,箱體壓力變化與外部激勵大小有直接關系。Kassemi 等[11]采用數值模擬分析橫向加速度引起的晃動對貯箱傳熱和壓降的影響。Montsarrat[12]以儲罐壓降為重點,預測了液氫貯箱中的晃動熱力學性能。李培昌等[13]研究發現貯箱內氣液界面在晃動激勵作用下被擾動,導致氣相冷凝增加,貯箱內壓力迅速降低等現象。Berglund 等[14]對液氫貯箱的相關預測表明,液氫可能被注入排氣口和泄放系統,導致推力不平衡和控制失效。2020年12月9日,Space X 對星艦運載器一枚原型試驗箭進行了首次高空飛行測試。但在著陸過程中,由于頭部燃料貯箱壓力過低,導致著陸速度過快,飛船發生爆炸[15]。因此,為保證低溫貯箱的安全運行,有必要對飛行器用貯箱內液氫晃動熱力學響應進行深入研究。

本文針對飛行器用160 L 高真空多層絕熱液氫貯箱,數值模擬研究不同的晃動激勵、初始液體溫度和初始充滿率下的晃動熱力學特性,旨在為飛行器用低溫液氫貯箱防晃設計及優化提供技術參考。

2 計算模型

2.1 VOF 模型

液體晃動過程中,液氫的流動是一個自由界面運動的問題,且液氫貯箱內始終存在氣液界面,本文采用流體體積方法(Volume of Fluid,VOF)模型捕捉對相界面進行追蹤。VOF 方法是一種在固定的歐拉網格下的表面追蹤方法,其根據相份額φq的值來確定自由界面,如式(1)所示:

界面質量、動量以及能量平衡通過加載到離散項中的源項實現,用于求解該過程的控制方程。連續性方程如式(2)所示:

式中:ρ為密度,ui為平均速度矢量,Sm為源項。

動量方程如式(3)所示:

式中,ρ為密度;P為靜壓;τij為應力張量;gi和Fi分別為i 方向上的重力體積力和外部體積力,Fi包含其他模型相關源項。應力張量由式(4)給出:

式中,τij為應力張量;μ為動力粘度;ui、uj分別為i、j 方向上的平均速度矢量。

能量方程如式(5)所示:

式中,ρ為密度;u為流體速度;T為溫度;cp為比熱容;k為流體換熱系數;gradT為溫度梯度,用來表示溫度的變化率;Sh為流體的內熱源及由于粘性作用流體機械能轉換為熱能的部分,簡稱粘性耗散項。

在外部漏熱作用下,氣液界面會發生傳熱傳質過程。相變模型采用Lee 模型,該模型假設界面相變發生在恒壓和準熱平衡狀態下,如式(6)所示:

式中:Sm為質量能量源項;rl、rv分別為液相、氣相的氣化潛熱;αl、αv為液相、氣相所占相份額;Tl、Tv為液相、氣相溫度;Tsat為飽和溫度。

2.2 數值模型

研究所用液氫貯箱為球形貯箱,主要由內容器和外殼組成。因內容器中液體晃動是本文重點研究內容,為簡化問題,僅針對內容器建模。圖1 展示了本文所研究液氫貯箱,其內容器直徑為700 mm,壁厚為2 mm,罐壁材料為5083 鋁合金。貯箱外部環境溫度為293 K,與低溫流體間存在較大溫差,部分環境漏熱q將滲入低溫貯箱。為研究貯箱內流體晃動問題,對貯箱施加水平晃動激勵,并在箱體內部氣相區設有壓力監測點以研究貯箱內流體動力特性,具體如圖1 所示。為節省計算資源,采用二維模型計算求解。

圖1 低溫液氫貯箱內罐示意圖Fig.1 Schematic diagram of the inner tank of a cryogenic liquid hydrogen tank

2.3 網格劃分

采用ANSYS 自帶前處理軟件Meshing 對液氫貯箱進行網格劃分。計算采用結構化網格,近壁處進行加密處理,網格模型及網格無關性驗證如圖2 所示。飛行器機動過程中晃動激勵為0.11 m/s 時,選擇了3 種不同數量(70 967、115 489 與173 892)的計算網格進行網格無關性驗證。對比不同計算網格下的晃動力變化可知,3種計算網格下箱體所受晃動力幾乎完全重合,晃動力波動趨勢存在微小差異,表明網格數量對箱體所受晃動力的影響較小。因此,為減少計算資源并提高計算精度,選擇數目為115 489 的網格進行后續數值模擬研究[16]。

圖2 網格模型及無關性驗證Fig.2 Grid model and grid independence verification

2.4 初始條件

Agui 和Arndt 等[17-18]針對低溫杜瓦內所做實驗和數值結果表明,貯箱內的實際溫度分布呈熱分層狀態,即氣液相間無溫度階躍。為保證貯箱發射前氣液相間溫度良好的過渡,液相初始溫度為20.0 K,在同一初始充滿率不同晃動激勵、不同初始液體溫度及不同初始充滿率下的氣枕區域從界面處(20.0～21.5 K)到貯箱頂部(25 K)采用線性溫度分布如圖3 所示。罐內初始充滿率50%,即液位高度 350 mm,初始氣枕壓力0.15 MPa。

圖3 液氫貯箱內初始溫度分布Fig.3 Initial temperature distribution in the liquid hydrogen tank

2.5 邊界條件

2.5.1 晃動激勵

為模擬飛行器用液氫貯箱晃動,選擇不同頻率不同振幅的簡諧激勵作為動量源項添加到計算模型中。相應的簡諧波動方程如式(7)所示:

式中,A為橫向激勵振幅,本文所取振幅分別為0.018 m、0.036 m、0.071 m;激勵頻率f=1.0 Hz;t為時間,計算步長為0.001 s。

貯箱與飛行器具有相同的飛行速度,可對貯箱壁面施加飛行器運行的速度邊界條件(后文簡稱晃動激勵),以模擬貯箱內部的液氫晃動,模擬設置中,在x軸方向(如圖1 所示)對貯箱罐壁施加如下速度晃動激勵,作用時間為2.5 s。

相應的速度晃動激勵如式(8)所示:

式中,v為速度激勵,通過位移激勵x求導得到,表示為x＇。

為實現液氫貯箱所受到的簡諧激勵,本文通過編寫用戶自定義函數(User Defined Function,UDF)并結合動網格技術對氣液界面波動過程進行有效預測。

2.5.2 漏熱熱流密度

外界熱流量進入貯箱的途徑,主要包括通過絕熱層、管路及玻璃鋼支撐結構的漏熱量。

通過絕熱層的漏熱量如式(9)所示:

式中,λMLI取0.000 22 W(m·K)為高真空多層絕熱表觀導熱系數,MLI 表示多層絕熱,A為絕熱層表面積,Ta取293 K,Td取20 K 分別為環境溫度及介質溫度,通過絕熱層漏熱量為3.66 W。

通過氣相管的漏熱量如式(10)所示:

式中,λ為5083 鋁無縫管熱導率,單位W/(m·K);Fe為氣相管的換熱面積,單位m2;Le為氣相管的管長,單位m;Ta、Td分別為環境溫度和介質溫度,單位K。

通過液相管的漏熱量如式(11)所示:

通過管路的綜合漏熱量為2.30 W,如式(12)所示:

通過玻璃鋼支撐結構的綜合漏熱量如式(13)所示:

式中,n為支撐結構數量;λs為玻璃鋼熱導率,單位W(m·K);A為玻璃鋼橫截面積,單位m2;L為玻璃鋼長度,單位m;Ta、Td分別為環境溫度和介質溫度,單位K。通過玻璃鋼支撐的綜合漏熱量為2.25 W。

通過絕熱層、管路和玻璃鋼支撐結構進入液氫貯箱的總漏熱量如式(14)所示:

式中,Q為液氫貯箱總漏熱量,8.21 W。通過內容器壁面的漏熱熱流密度:

式中,為內容器外表面的面積,1.56 m2。

數值模擬中,對漏熱熱流做時均化處理,均視為第二類邊界條件,其值為5.26 W/m2。

2.6 求解設置

采用Fluent 雙精度求解器對液氫儲罐內液氫晃動的熱力耦合特性進行瞬態數值模擬。貯箱內氣枕空間增壓氣體為氫氣,計算工質為液氫。氫氣采用理想氣體模型,液氫密度采用Boussinesq假設,運動粘度等參數均參考物性軟件NIST。參量的離散格式設置如表1 所示。

表1 離散格式Table 1 Discretization scheme

2.7 數值模型驗證

由于液氮在物理性質上與液氫最為接近,而流體的晃動與運動粘度直接相關,液氮與液氫的運動粘度較為接近,標況下二者比值如式(16)所示:

式中,υLN2為液氮運動黏度,υLH2為液氫運動黏度。

采用軟件ANSYS Fluent 19.0 模擬Lacapere[19]和Agui 等[17]針對立式低溫恒溫器所進行的液氮晃動實驗,以驗證本文所構建數值模型的有效性和準確性。模擬過程中的初始條件和邊界條件與實驗保持一致,計算處理與設置與本文所構建模型相一致。獲得了低溫恒溫器內氣枕壓力值。因液氮晃動實驗過程重點關注低溫恒溫器內部壓力變化,此處選取實驗值與模擬結果進行對比,兩者之間的相對誤差如圖4 所示。

圖4 實驗值與模擬值間的相對誤差Fig.4 Relative error between the experimental value and the simulated value

由圖4 可知,氣枕壓力降幅高于實驗,原因在于數值模擬未考慮氣液界面附近的液相熱分層。在40 s 的數值計算過程中,除前9 s 和后3 s 的模擬值與實驗值存在較大偏差(最大為18%)外,其余時間兩者誤差均小于5%。結果表明,本文所構建的數值模型可用于運載飛行器用液氫貯箱內液氫晃動熱力學響應的研究。

3 結果與討論

3.1 不同晃動激勵下晃動熱力性能

為研究不同晃動激勵對飛行器用液氫貯箱內熱力學性能的影響,比較了晃動激勵為0.00 m/s,0.11 m/s,0.22 m/s,0.44 m/s 時液氫貯箱內氣枕壓力的變化情況,氣枕壓力變化曲線及貯箱溫度場分布分別如圖5 和圖6 所示。

圖5 氣相測點Pv 處壓力變化曲線Fig.5 Pressure change curve at the gas phase measurement point Pv

由圖5 可知,當施加不同晃動激勵時,貯箱氣枕壓力并未立即迅速降低,隨晃動激勵的擾動,破壞了貯箱內熱力學平衡。①當晃動激勵為0 m/s時,因環境熱量漏入罐中,氣枕壓力出現先降后升的趨勢,最大壓降2 Pa,最大壓升15.4 Pa。②當晃動激勵為0.11 m/s 時,在2.5 s 的數值模擬過程中,氣液界面附近的兩相流只出現了小幅晃動,壓力曲線呈現出小振幅的線性響應,最大壓降120.3 Pa。③當晃動激勵為0.22 m/s 時,貯箱內液體表現出較強的非線性響應。結合溫度云圖6(c)可知,1.6 s 時液體以剪切流的形式到達氣枕空間,同時由于外界晃動激勵的作用,貯箱內流體往復運動,致使剪切流發生波浪破碎,產生大量飛濺液滴,氣枕壓力出現大幅下降的拐點,直至壓降速率達到最大后,貯箱內被擾動的熱力學狀態重新達到平衡,最大壓降為6084.5 Pa。④當晃動激勵為0.44 m/s 時,壓力曲線表現出高振幅的非線性響應。前1.1 s 氣枕壓力線性降低,1.3 s 后,液體波浪到達氣枕空間,并夾帶大量過熱氣體。相比0.22 m/s 而言,液體到達氣枕空間后產生更多量的飛濺液滴,導致氣、液相間接觸面積急劇增大,在增強氣、液兩相間換熱的同時伴隨著相變傳質,導致氣枕壓力大幅降低,氣枕空間最大壓降9158.3 Pa。

圖6 不同晃動激勵下溫度場分布云圖Fig.6 Temperature contour under different sloshing excitations

總之,氣枕壓力隨晃動的進行持續降低,且晃動激勵越大氣枕壓力降幅越大,晃動激勵為0.11、0.22、0.44 m/s 時的氣枕壓降分別為120.3、6084.5、9158.3 Pa。

由圖7 可知不同晃動激勵下氣液界面處測點Ti處溫度隨時間的變化規律。在未施加晃動激勵時,測點溫度微幅波動升高,最大增幅為0.05 K;晃動激勵為0.11 m/s 時,測點溫度小幅波動升高,最大增幅為0.315 K;晃動激勵為0.22 m/s 時,在2.08 s 之前,測點溫度波動升高,最大增幅為1.809 K,在2.08 s 后,受流動粘性及阻尼運動的影響,測點溫度隨時間波動降低。晃動激勵為0.44 m/s 時,與晃動激勵為0.22 m/s時測點溫度的波動變化趨勢相似,在1.91 s 之前,測點溫度波動升高,最大增幅為1.392 K;隨后受流體粘性及阻尼運動的影響,測點處溫度波幅逐漸減小。

圖7 不同晃動激勵下氣液界面測點Ti 溫度變化Fig.7 Temperature change at point Ti of gas-liquid interface under different sloshing excitations

由以上分析可知,在2.5 s 的數值計算時間內,不同晃動激勵下氣液界面處測點Ti處的溫度隨時間波動增大,且晃動激勵越大,波動越劇烈。原因在于當施加不同晃動激勵時,貯箱軸對稱上的測點受氣相傳熱及壁面自然對流作用,被加熱液體會向氣液界面處流動積聚,因此,氣液界面處測點溫度隨時間波動升高。

3.2 不同初始液體溫度下氣枕壓力變化

當初始氣枕壓力為150 kPa 時,液氫的飽和溫度為21.77 K,選擇3 種初始液體溫度分別為20.0、21.0、21.5 K,其值均小于液氫的飽和溫度,即該溫度下的液相皆處于過冷狀態。

晃動激勵為0.22 m/s、初始充滿率為50%和初始液體溫度分別為20.0、21.0、21.5 K 時氣相測點Pv的氣枕壓力變化如圖8 所示。

圖8 不同初始液體溫度下氣相測點Pv 處壓力變化Fig.8 Pressure changes at the gas phase measurement point Pv under different initial liquid temperatures

由圖8 可知,在過冷液體冷凝下,前1.5 s 不同初始液體溫度下的氣枕壓力隨時間線性降低,1.5 s 后剪切流到達氣枕空間,增大了氣液相間接觸面積。此外,對液氫而言,初始液體溫度越低,過冷度越大,當晃動激勵作用于罐壁,使得氣液相間接觸面積增大時,氣枕空間過熱氣相被初始液體溫度較低的液氫極大冷凝,致使氣枕空間發生更大壓降。因此,氣枕壓降隨初始液體溫度的降低而增大。當初始液體溫度從20.0 K 升高到21.5 K 時,測點處的最終氣枕壓力值分別為143 921、144 752、146 098 kPa,相對應的壓降分別為6079、5248 和3902 Pa。

3.3 不同初始充滿率下氣枕壓力變化

為探究晃動激勵為0.22 m/s、初始液體溫度為20.0 K 和初始充滿率分別為30%、40%、50%、60%、70%下的晃動熱力性能,在罐體對稱軸上設置了氣相測點Pv(0,600)。該測點在不同初始充滿率、同一晃動激勵下的氣枕壓力變化曲線如圖9 所示。

圖9 不同初始充滿率下氣相測點Pv 處壓力變化Fig.9 Pressure changes at the gas phase measurement point Pv under different initial full rates

由圖9 可知,不同初始充滿率下氣枕壓力不斷降低,且初始充滿率越高,過冷液體中貯存了更多的冷量,氣枕壓力降幅越大。隨著初始充滿率從30%增加到70%,氣枕壓力從150 kPa 下降到148 095、146 242、143 915、143 523.83 和141 661 kPa,相應的壓降分別為1905、3758、6085、6476 和8339 kPa。

4 結論

采用CFD 技術數值模擬研究了外部晃動激勵下飛行器用液氫貯箱內低溫推進劑晃動熱力學響應分析,結論如下:

1)晃動增強了過熱氫氣對過冷液體和罐壁的傳熱,晃動激勵越大,氣枕空間壓降越大。

2)初始液氫溫度對氣枕壓力有較大影響,氣枕壓降隨初始液體溫度的降低而增大。

3)初始充滿率越高,過冷液體中貯存的冷量越大,氣枕壓力對界面的冷凝更敏感,氣枕壓降越大。

4)貯箱內壓降受晃動激勵、初始液氫溫度及初始充滿率的影響,其中外部晃動激勵極大的擾動了液氫貯箱內氣液界面處的熱力學平衡,導致貯箱氣枕壓力大幅降低,影響飛行器的穩定運行。因此,亟需采取合理的增壓或防晃措施來維持貯箱壓力的穩定性,為飛行器的安全運行提供保障。