微重力環(huán)境下低溫推進(jìn)劑貯箱內(nèi)氣液界面形變特性研究

閆春杰，鄭永煜，楊 祺，楊 鵬，劉迎文，王小軍

（1.蘭州空間技術(shù)物理研究所 真空技術(shù)與物理重點實驗室，蘭州 730000；2.西安交通大學(xué) 能源與動力工程學(xué)院熱流科學(xué)與工程教育部重點實驗室，西安 710049）

0 引言

低溫液體推進(jìn)劑在加注和航天器飛行過程中分別受到壓降和加熱等的影響，以致部分蒸發(fā)，在貯箱內(nèi)部形成“氣枕”區(qū)?！皻庹怼眳^(qū)內(nèi)比熱容較低的氣體吸收外界熱量造成貯箱內(nèi)的壓力急劇升高并產(chǎn)生波動，威脅低溫液體存儲系統(tǒng)的安全運行。通常采用熱力學(xué)排氣、主動制冷與被動絕熱相結(jié)合的零蒸發(fā)存儲方法控制低溫貯箱內(nèi)的壓力，降低低溫液體的蒸發(fā)量。對于常規(guī)地面工況，氣液的空間構(gòu)型清晰，容易實現(xiàn)直接排氣的壓力控制策略，不會出現(xiàn)不必要的低溫液體無效損失。但對處于空間微重力環(huán)境下的低溫推進(jìn)劑，氣液空間構(gòu)型未知，在“氣枕”形狀和位置不確定時，簡單的排氣放壓方式容易造成液體的排放損失[1]。此外，微重力條件下氣液空間構(gòu)型嚴(yán)重地影響主動控壓技術(shù)的控壓效率。因此，為保證低溫推進(jìn)劑在軌安全儲存，提高主動控壓技術(shù)的控壓效率，開展微重力下“氣枕”的構(gòu)型及形變特性的研究至關(guān)重要。

國內(nèi)外研究人員圍繞低溫貯箱內(nèi)“氣枕”形狀開展了研究。Charles等[2-3]、Mohammad等[4]將低溫貯箱內(nèi)的“氣枕”做了集中處理，仿真結(jié)果表明，受到殘余重力的影響，“氣枕”區(qū)域沿著與重力加速度相反的方向從貯箱中部逐步升至貯箱頂部位置，整個過程歷時567 s；隨后，他們又模擬了50%、90%兩種充灌率下“氣枕”由中心升至頂部及到達(dá)頂部后75 d時間內(nèi)貯箱內(nèi)低溫液體溫度場、流場和壓力變化，并利用零蒸發(fā)儲存技術(shù)成功控制貯箱內(nèi)部壓力。Thornton等[5-6]圍繞氣液界面上能量的進(jìn)出、界面的變形和噴射液體對自由液面形狀的改變等問題對箱內(nèi)液體流動狀況的影響進(jìn)行了深入研究，仿真結(jié)果表明，噴射的液體會將自由液面從初始的橢球分割成兩半。李章國等[7]利用連續(xù)表面張力模型討論了在軌航天器穩(wěn)定運行狀態(tài)下低溫液氨在貯箱內(nèi)的形態(tài)分布，分析討論了重力加速度、液-固接觸角等因素對低溫液氨氣液界面構(gòu)型的影響，仿真結(jié)果顯示，邦德數(shù)Bo較大時，重力對氣液界面構(gòu)型起主導(dǎo)作用；Bo較小時，液-固接觸角決定氣液界面構(gòu)型。鄧新宇等[8]利用Flow-3D軟件驗證推進(jìn)劑管理過程中的重定位問題，研究了加速度、推進(jìn)劑充灌率初始形態(tài)等因素對重定位過程的影響，計算結(jié)果表明，推進(jìn)劑全部位于貯箱前底時重定位的工況最惡劣，在相同沉底推力下，充灌率越高重定位時間越短。魏延明等[9]研究了微重力環(huán)境下低溫貯箱內(nèi)推進(jìn)劑氣液界面的靜平衡和重定位問題，利用相似理論縮比模型進(jìn)行了一系列的數(shù)值計算，并與落塔實驗對比，測算出重定位的時間，為表面張力貯箱的設(shè)計提供了依據(jù)。肖立明等[10]運用VOF法對微重力下液氧貯箱內(nèi)自由液面變形問題進(jìn)行模擬，仿真得到的液體爬升高度和表面張力波傳播速度與理論計算結(jié)果的符合性良好，驗證了VOF方法的可行性。劉趙淼等[11]同樣采用VOF法研究了不同重力加速度和液-固接觸角對板式表面張力貯箱內(nèi)液體流動過程的影響。張鎧等[12]以液氫為對象，研究了不同重力加速度、液-固接觸角和充灌率等因素對貯箱內(nèi)氣液界面構(gòu)型的影響，計算結(jié)果表明，液-固接觸角對氣液界面構(gòu)型的影響占主導(dǎo)作用，充灌率的影響不大。

國外對微重力情況下低溫推進(jìn)劑儲存系統(tǒng)“氣枕”形變特性的研究情況極為保密，尚未查閱到可供參考的系統(tǒng)的研究數(shù)據(jù)。著眼于未來空間應(yīng)用對低溫液體儲存技術(shù)的需求，亟需開展低溫推進(jìn)劑貯箱內(nèi)“氣枕”形變特性的研究。本文采用VOF法結(jié)合連續(xù)表面張力（CSF）模型[10-14]，研究貯箱內(nèi)“氣枕”構(gòu)型以及重力水平、推進(jìn)劑液-固接觸角等因素對于氣液界面構(gòu)型的影響，為空間推進(jìn)劑存儲系統(tǒng)的低溫流體管理提供理論依據(jù)。

1 數(shù)學(xué)物理模型

1.1 物理模型

低溫液體貯箱的結(jié)構(gòu)示意圖如圖1所示，貯箱由中間的圓柱段和兩端的橢圓封頭構(gòu)成，其外壁面鋪設(shè)有多層絕熱材料，以減少外界進(jìn)入貯箱的熱量。貯箱內(nèi)部包含導(dǎo)熱元件和噴射泵，導(dǎo)熱元件可將熱量從壁面導(dǎo)入貯箱中心區(qū)域，噴射泵起攪混貯箱內(nèi)流體的作用，兩者協(xié)同工作能實現(xiàn)控制貯箱壓力以及消除熱分層現(xiàn)象的目的。采用的低溫液體為液氧，貯箱體積為3.6 m3，圓柱段外徑A為1.4 m，高度B為2.6 m；導(dǎo)熱元件高度C為1.3 m；噴射泵高度D為0.55 m；橢圓封頭短軸長E為0.75 m，長軸長F為1.4 m。

圖1 液氧低溫貯箱結(jié)構(gòu)示意圖Fig.1 Schematic diagram of liquid oxygen cryogenic storage tank

由于本文研究重點為“氣枕”的形變特性以及影響氣液界面構(gòu)型的主要因素，因此在建模過程針對箱內(nèi)的布置做出簡化：（1）不考慮噴射泵以及導(dǎo)熱元件對界面構(gòu)型的影響；（2）由于界面構(gòu)型時間遠(yuǎn)小于熱量從壁面導(dǎo)入低溫液體內(nèi)部所需時間，因此不考慮冷卻和外界熱負(fù)荷的影響；（3）模擬“氣枕”構(gòu)型時，初始時刻將“氣枕”區(qū)域設(shè)置為平坦的平面。

1.2 數(shù)學(xué)模型

采用VOF法追蹤微重力下貯箱內(nèi)“氣枕”的位置變化，通過對整個流體區(qū)域求解質(zhì)量守恒、動量守恒與能量守恒控制方程，確定液氧在貯箱內(nèi)的界面構(gòu)型。

VOF模型用流體體積函數(shù)f來表示相函數(shù)，流體體積函數(shù)f的定義為每個空間單元內(nèi)，流體所占據(jù)的體積與該單元最大可容納的體積之比。流體體積函數(shù)f的控制方程如式（1）：

液體與固體的液-固接觸角θ反映液體分子與固體分子間的附著力與液體分子之間內(nèi)聚力的相對大小，表征液體對固體面的浸潤性。液-固接觸角滿足經(jīng)典的Yong氏方程：

式中：σgs為氣-固界面張力；σgl為氣-液界面張力；σls為液-固界面張力。

在微重力環(huán)境下，由于表面張力作用的權(quán)重增加，氣-液相界面發(fā)生彎曲，界面曲率及壓力差滿足Laplace-Young表面張力方程：

式中：Δp為彎曲界面內(nèi)外壓力差，Pa；R1與R2為自由液面任一點上相互垂直的正截面的曲率半徑，m。

微重力環(huán)境下，流體動力學(xué)特性與常重力下存在顯著區(qū)別。隨著重力因素的逐漸消退，表面張力成為影響氣液界面（或液體）形狀、位置的主要因素。衡量表面張力（毛細(xì)作用）與重力因素相對大小的準(zhǔn)則數(shù)為邦德數(shù)Bo：

式中：R為低溫貯箱的特征尺寸，m。

2 數(shù)值計算方法

選擇適用于瞬態(tài)流動的PISO方案處理速度與壓力的耦合，選擇PRESTO算法計算壓力插值，選擇幾何重建（Geo-Reconstruct）格式捕捉自由界面附近的插值，選擇增強壁面函數(shù)處理低雷諾數(shù)下的復(fù)雜回轉(zhuǎn)運動，選擇壁面無穿透、無滑移邊界條件。

計算過程中，設(shè)定氣體為理想氣體，液體為不可壓縮流體。在VOF法CSF模型中，為了提高解的穩(wěn)定性，在多相模型中設(shè)置可壓縮的氣相為主相，液相為次相，用于平衡可壓縮相的壓力梯度和表面張力之間的相互作用，加速整個計算的收斂過程。因此，當(dāng)網(wǎng)格單元中體積函數(shù)f=1時，該單元為主相，即網(wǎng)格單元全部被氣相所占據(jù)；當(dāng)f=0時，該單元全部為液相；f介于0～1之間表明該網(wǎng)格處為自由液面。

有資料[15]表明，將低溫貯箱發(fā)射升空至界面穩(wěn)定階段全部耗時約500 s，而熱量由壁面導(dǎo)入液氧所需的時間量級為td=Rd2/al=6.3×106s，式中Rd為貯箱的當(dāng)量半徑，m；al為熱擴(kuò)散系數(shù)，m2·s-1，其數(shù)值遠(yuǎn)大于500 s。因此，可認(rèn)為熱量的擴(kuò)散停留在壁面附近很小的區(qū)域內(nèi)，基本不會影響到內(nèi)部流體的溫度分布以及“氣枕”形狀的動態(tài)變化。在“氣枕”穩(wěn)定之前貯箱內(nèi)低溫液體的溫升約為（以充灌率95%為例）DT=Q/mcp=5.46×10-4K，式中Q為進(jìn)入貯箱的熱量，J；m為貯箱內(nèi)流體質(zhì)量，kg；cp為流體的比熱，kJ·kg-1·K-1，如此微小的溫升可認(rèn)為在“氣枕”形變過程中，工質(zhì)物性例如表面張力、黏度等熱物性參數(shù)基本保持不變。因此，在建模過程中，將其設(shè)置為定值，如表1所列。

表1 物性設(shè)置Tab.1 Physical property setting

為保證計算的準(zhǔn)確性，設(shè)置計算時間步長為0.001 s，計算過程中取較小的松弛因子以確保計算收斂。收斂判據(jù)為：第一，某次迭代的相對誤差小于預(yù)設(shè)的容許誤差，其中連續(xù)性殘差小于10-3，其余殘差小于10-6；第二，在某次迭代中指定點的溫度值或者速度值不再隨時間發(fā)生變化。

3 計算結(jié)果與討論

氣液界面構(gòu)型取決于重力水平（表面張力）、液-固接觸角以及充灌率等因素。從文獻(xiàn)[12]的結(jié)論可知，充灌率對氣液界面構(gòu)型影響不大，因此本文重點研究液氧貯箱內(nèi)“氣枕”形變過程以及不同貯箱重力水平和液-固接觸角對其構(gòu)型的影響。

3.1 “氣枕”構(gòu)型

圖2為10-5g重力、5°液-固接觸角以及95%充灌率條件下液氧貯箱內(nèi)氣液自由界面隨時間的變化過程。其中，紅色部分表示氣相區(qū)域，藍(lán)色部分表示液相區(qū)域，貯箱內(nèi)Bo=5.37。由圖可知，在0～100 s過程中，液體的表面勢能轉(zhuǎn)化成位勢能，液體沿壁面依附爬升，氣液界面曲率半徑不斷減小。100 s時，氣液界面從初始的平面逐漸轉(zhuǎn)變成類似于橢球體的泡狀回轉(zhuǎn)體，液體完全覆蓋在貯箱的內(nèi)壁面上。在200～300 s過程中，氣液界面在液體內(nèi)部慣性的作用下發(fā)生輕微波動。300 s時，“氣枕”形成固定的回轉(zhuǎn)體型，如圖2（h）所示。因此，在微重力環(huán)境中，液體對貯箱壁壁面浸潤性良好的條件下，表面張力將驅(qū)動液體沿著壁面依附爬升直至最終完全浸潤壁面，而氣液界面也隨著液體的爬升形成穩(wěn)定的橢球狀回轉(zhuǎn)體。

圖2 10-5g重力、5°液-固接觸角、95%充灌率時液氧貯箱內(nèi)液面變化過程Fig.2 Liquid oxygen storage tank liquid level change process under 10-5g gravity，5°contact angle，95% filling rate

為了更清晰地反映“氣枕”動態(tài)形成過程中貯箱內(nèi)液體的流動狀態(tài)與速度分布，選取兩個特殊位置點監(jiān)測液體流動速度。其中，監(jiān)測點1設(shè)置在“氣枕”成型穩(wěn)定后正下方中心點，監(jiān)測點2設(shè)置在橢球封頂與柱形筒壁交匯處，兩個監(jiān)測點測到的“氣枕”形變過程中液體的速度變化如圖3所示。由圖可知，100 s時，監(jiān)測點1和2液體的流動速度接近零，說明液體表面勢能完全轉(zhuǎn)換為位勢能，液體完全覆蓋貯箱壁面。此后，由于低溫液體爬升運動的慣性，自由液面發(fā)生一定的波動，并且由于黏性耗散的存在，速度動能逐漸耗散成內(nèi)能。300 s后，兩監(jiān)測點液體流動速度均接近0，“氣枕”構(gòu)型達(dá)到穩(wěn)定，形成類似于橢球形的氣液分界面，結(jié)果與圖2一致。

圖3 “氣枕”形變過程兩監(jiān)測點液體流動速度變化曲線Fig.3 Velocity changes at two monitoring points in the process of gas-liquid interface deformation

3.2 液-固接觸角對“氣枕”形狀的影響

為了定量研究液-固接觸角對“氣枕”形變特性的影響，以95%充灌率的液氧貯箱為例，研究在10-5g重力環(huán)境下，液-固接觸角變化所引起的低溫液體貯箱內(nèi)氣液界面構(gòu)型的變化。其中，液-固接觸角變化范圍為0°～180°，由流體物性（如表面張力）和壁面的物理性質(zhì)（例如材質(zhì)和表面粗糙度）共同決定。液-固接觸角越小，液體潤濕性能越好。當(dāng)液-固接觸角為0°時，液體完全潤濕固體表面；當(dāng)液-固接觸角為180°時，液體完全不潤濕固體表面。為了便于分析，引入液體浸潤度Rin：

式中：Sin為“氣枕”穩(wěn)定時液體浸潤貯箱壁面的面積，m2；Sw為貯箱內(nèi)表面積，m2。液體浸潤度為1表明液體完全浸潤壁面，用浸潤度值可以確定推進(jìn)劑對貯箱壁面的覆蓋狀況，反映“氣枕”的構(gòu)型變化。

圖4為10-5g重力環(huán)境下，充灌率為95%時，貯箱內(nèi)液氧的浸潤度隨液-固接觸角的變化規(guī)律。由圖可知，當(dāng)液-固接觸角為0時，浸潤度為1，液氧完全浸潤壁面；當(dāng)液-固接觸角小于20°時可以實現(xiàn)0.98以上的良好浸潤性。隨著液-固接觸角不斷增大，液氧與壁面的附著力減弱，液體浸潤性不斷下降。當(dāng)液-固接觸角為150°時，浸潤度為0.861，小于平坦表面的浸潤度0.874，表明此時氣液界面呈現(xiàn)凸形，浸潤面積比液氧面水平狀態(tài)下（常重力g=9.8 m/s2工況）小。由于熱量直接通過“氣枕”進(jìn)入貯箱造成的壓力上升十分劇烈，因此應(yīng)將液-固接觸角控制在20°范圍內(nèi)，使液體的浸潤度大于0.98，以降低壓力上升速率，延長低溫液體在軌儲存時間。

圖4 10-5g重力下，充灌率95%時，貯箱內(nèi)液氧浸潤度變化Fig.4 Wettability change under the 10-5g gravity，95% filling rate in the liquid oxygen tank

3.3 重力水平對“氣枕”形狀的影響

充灌率為95%、液-固接觸角為5°時，貯箱內(nèi)液氧浸潤度隨重力水平變化規(guī)律如圖5所示。

由圖5可知，重力水平由10-5g增大至1g的過程中，液體浸潤度呈現(xiàn)下降趨勢。當(dāng)重力水平小于10-3g時，液體浸潤度隨著重力水平的下降線性上升，在此區(qū)間內(nèi)液體的爬升能力隨重力水平的下降迅速提升，表明表面張力成為決定氣液界面構(gòu)型的主要因素。當(dāng)重力水平大于10-3g時，液體爬升能力較弱，液體浸潤度隨重力水平的上升呈現(xiàn)出緩慢下降的趨勢，氣液界面的構(gòu)型逐漸表現(xiàn)出類似于地面的情形，說明在此區(qū)間內(nèi)重力水平成為決定氣液界面構(gòu)型的主要因素。通常情況下，近地軌道的重力水平在10-6g～10-4g區(qū)間，由上述模擬結(jié)果可知，“氣枕”的構(gòu)型為類似于橢球形的彎曲界面，液體能較好地浸潤貯箱內(nèi)壁，有利于貯箱內(nèi)低溫推進(jìn)劑的存儲及壓力的有效控制。

圖5 充灌率95%、液-固接觸角5°時貯箱內(nèi)液氧浸潤度變化Fig.5 Liquid oxygen wettability in the tank when the filling rate is 95% and the liquid-solid contact angle is 5°

4 結(jié)論

本文采用VOF法結(jié)合CSF模型，研究貯箱內(nèi)“氣枕”構(gòu)型以及重力水平、液-固接觸角等因素對液氧貯箱內(nèi)氣液界面構(gòu)型的影響，得到了如下結(jié)論：

（1）在微重力環(huán)境中，由于表面張力作用的存在，液體沿著壁面依附爬升，貯箱內(nèi)的“氣枕”最終構(gòu)型為橢球形回轉(zhuǎn)體。

（2）在10-5g重力水平下，貯箱內(nèi)液體的浸潤性隨液-固接觸角增大而減小。為了滿足0.98以上的液體浸潤度，在一定表面張力下，應(yīng)使用液-固接觸角小于20°的貯箱材料，以增大浸潤面積，延長低溫推進(jìn)劑的在軌儲存時間。

（3）當(dāng)重力水平位于10-4g～10-3g區(qū)間內(nèi)時，液體的爬升能力迅速提升，重力因素與表面張力的權(quán)重在此區(qū)間內(nèi)發(fā)生顯著變化。近地軌道的重力水平在10-6g～10-4g區(qū)間，貯箱內(nèi)“氣枕”的構(gòu)型為類似于橢球形的彎曲界面，液體能較好地浸潤貯箱內(nèi)壁，有利于低溫推進(jìn)劑的長期存儲。