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

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

1 引言

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

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

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

2 計算模型

2.1 VOF 模型

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

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

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

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

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

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

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

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

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

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

2.2 數(shù)值模型

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

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

2.3 網(wǎng)格劃分

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

圖2 網(wǎng)格模型及無關(guān)性驗證Fig.2 Grid model and grid independence verification

2.4 初始條件

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

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

2.5 邊界條件

2.5.1 晃動激勵

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

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

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

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

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

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

2.5.2 漏熱熱流密度

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

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

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

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

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

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

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

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

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

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

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

式中,為內(nèi)容器外表面的面積,1.56 m2。

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

2.6 求解設(shè)置

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

表1 離散格式Table 1 Discretization scheme

2.7 數(shù)值模型驗證

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

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

由圖7 可知不同晃動激勵下氣液界面處測點Ti處溫度隨時間的變化規(guī)律。在未施加晃動激勵時,測點溫度微幅波動升高,最大增幅為0.05 K;晃動激勵為0.11 m/s 時,測點溫度小幅波動升高,最大增幅為0.315 K;晃動激勵為0.22 m/s 時,在2.08 s 之前,測點溫度波動升高,最大增幅為1.809 K,在2.08 s 后,受流動粘性及阻尼運動的影響,測點溫度隨時間波動降低?；蝿蛹顬?.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 的數(shù)值計算時間內(nèi),不同晃動激勵下氣液界面處測點Ti處的溫度隨時間波動增大,且晃動激勵越大,波動越劇烈。原因在于當(dāng)施加不同晃動激勵時,貯箱軸對稱上的測點受氣相傳熱及壁面自然對流作用,被加熱液體會向氣液界面處流動積聚,因此,氣液界面處測點溫度隨時間波動升高。

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

當(dāng)初始?xì)庹韷毫?50 kPa 時,液氫的飽和溫度為21.77 K,選擇3 種初始液體溫度分別為20.0、21.0、21.5 K,其值均小于液氫的飽和溫度,即該溫度下的液相皆處于過冷狀態(tài)。

晃動激勵為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 后剪切流到達(dá)氣枕空間,增大了氣液相間接觸面積。此外,對液氫而言,初始液體溫度越低,過冷度越大,當(dāng)晃動激勵作用于罐壁,使得氣液相間接觸面積增大時,氣枕空間過熱氣相被初始液體溫度較低的液氫極大冷凝,致使氣枕空間發(fā)生更大壓降。因此,氣枕壓降隨初始液體溫度的降低而增大。當(dāng)初始液體溫度從20.0 K 升高到21.5 K 時,測點處的最終氣枕壓力值分別為143 921、144 752、146 098 kPa,相對應(yīng)的壓降分別為6079、5248 和3902 Pa。

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

為探究晃動激勵為0.22 m/s、初始液體溫度為20.0 K 和初始充滿率分別為30%、40%、50%、60%、70%下的晃動熱力性能,在罐體對稱軸上設(shè)置了氣相測點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,相應(yīng)的壓降分別為1905、3758、6085、6476 和8339 kPa。

4 結(jié)論

采用CFD 技術(shù)數(shù)值模擬研究了外部晃動激勵下飛行器用液氫貯箱內(nèi)低溫推進劑晃動熱力學(xué)響應(yīng)分析,結(jié)論如下:

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

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

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

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