重力跳變引起的貯箱內氣液界面波傳播規律*

魏 列 杜王芳 薛子揚 何發龍 ,4 李 凱 趙建福

1(中國科學院力學研究所 微重力重點實驗室 北京 100190)

2(中國科學院大學工程科學學院 北京 100049)

3(北京工業大學材料與制造學部 北京 100124)

4(天津城建大學能源與安全工程學院 天津 300384)

0 引 言

隨著航天技術的發展，航天任務的復雜性不斷提高，對液體火箭發動機多次關機重啟的要求日益普遍，相應對推進劑管理系統的要求也越來越嚴格[1-3]。為了確保液體火箭發動機在關機滑行段之后再次啟動操作能夠順利進行，必須考慮貯箱內氣液界面如何響應重力或加速度水平改變的問題。由于推進劑在氣液兩相狀態下密度差異巨大，不同重力條件下氣液兩相分布特征迥異，重力變動時氣液界面和氣液兩相介質的運動一直都是先進空間流體管理技術研發的關鍵[4]。在液體火箭發動機關機滑行段航天器處于微重力條件下，體積力會減弱甚至完全被抑制，這導致控制流體行為的主要因素變為毛細力，液態推進劑不能像在重力場中那樣自發沉降到貯箱底部。液體火箭發動機再次啟動之前，往往需要使用輔助推力器將推進劑重新定位到貯箱供液口處。因此，研究推進劑在重力變動之后的界面動力學行為具有直接的應用價值[5]。

Bond（Bo）數是影響不同重力條件下貯箱內氣液兩相分布及相應流動特征的關鍵參數之一[6,7]。對于完全潤濕的液體（靜態接觸角為0°），直立圓管容器中穩定的氣液界面所對應的臨界Bo數（重力變化不影響容器中彎月面形狀穩定性的最大數值）是0.84[6-9]。Reynolds和Satterlee[10]給出了球形容器中彎月面的臨界Bo數的解析表達式。Weislogel和Hsieh[11]給出了直立方柱體內臨界Bo數隨接觸角的變化，接觸角為90°時數值最大。Bretherton[8]采用潤滑流假設求得0.84＜Bo＜1.05時氣泡中心點的遷移速度。Davies和Taylor[12]利用無黏勢流理論求出Bo＞＞1時倒置直立圓管中氣泡的運動速度。Masica和Petrasb[5]利用低重力和常重力實驗數據擬合出了Bo＞1時豎直圓管中氣泡遷移速度的經驗關聯式，其中忽略了黏性對氣泡運動的影響。

為了考慮黏性力的影響，可以利用Bo數和Reynolds數（Re）來表示直立圓管中氣泡的遷移速度[13]。Labus和Masica[14]定性地描述了球形貯箱中兩種初始位形的自由界面在重定位過程中的動態演變特征。Masica和Salzman[15]利用縮比的半人馬貯箱研究了初始彎曲的自由界面在重定位過程中的涌泉現象，發現區分界面涌泉是否發生的臨界Weber數（We=ρU2R/σ）是4。Salzman等[16]用落塔實驗研究了縮比半人馬液氫貯箱內液體的重定位過程，并估算了第一次收集貯箱頂部液體所需的時間。由于界面涌泉會在液體匯聚時形成，因而難以預測完整的沉底過程所需的時間。Hochstein等[17]數值模擬了幾何相似的貯箱在相同充液率下的重定位過程，提出了一個新的無量綱參數Settling數（Se=μ(Rg)1/2/σ）來表征黏性力、重力與表面張力間的相對作用強度，得到貯箱在重定位的過程中的速度增量正比于Se數。為了減少沉底推力所使用的推進劑，Patag等[18]采用了脈沖式沉底推力的方法，并定義了一個脈沖式沉底過程中的無量綱時間。通過數值模擬的方式發現，當無量綱時間小于并接近1.4時，其沉底操作的效果是最優的。Hochstein等[19]通過假定液體在重定位過程中作剛體運動，導出了一個能夠有效預測不同充液率或尺寸的貯箱完成重定位時航天器速度增量的無量綱參數，但只適用于Bo＞10的情況。

上述研究關注于貯箱內施加沉底推力時氣液界面的運動及氣液兩相流動特征，即重力由無向有變化的過程，而較少關注突然失重如關機過程中的演變。在后一情形中，由于重力突然減弱或消失，毛細力會驅動接觸線運動，進而引起界面的波動，并演變到新的平衡態穩定下來。直立圓柱容器中液體的平衡位形取決于Bo數和液體接觸角[20]。但在復雜幾何中無法只利用Bo數和接觸角來預測容器中的界面位形，需要采用系統表面能最小原理來進行預測[21]。Siegart等[22]進一步給出了球形、圓柱和環形容器內氣液界面到達平衡界面所需時間的關聯式，驗證了We數能有效預估貯箱內自由界面到達平衡位形所需的時間。在自由界面重新平衡之前，貯箱內的液面會經歷一個阻尼振蕩的過程，而低重力下氣液界面的共振頻率和阻尼率取決于接觸線邊界條件[23]。Kaukler[24]實驗研究了潤濕轉變溫度對重力突降過程中自由界面振蕩的影響，發現接觸線摩擦使得部分潤濕液體的周期在振蕩過程中不斷增加。Dreyer等[25-28]利用數值和實驗研究了重力水平變小后自由界面的運動及界面振蕩的阻尼特性。Li等[29]用落塔實驗和數值模擬研究了微重力條件下水作為工質的自由液面的振蕩頻率與貯箱尺寸和充液比的關系。Li等[30,31]采用動態接觸角模型數值研究了貯箱內自由界面在不同充液比以及不同重力水平下的振蕩特性。Xiao等[4]數值研究了重力突降后貯箱內界面波動演化特征。在此基礎上，Wei等[32]進一步研究并驗證了原型和縮比模型之間的運動相似性，但未能系統地研究和分析不同Bo數下貯箱內各作用力之間的關系。

上述研究有助于理解貯箱內流體運動特性，但對貯箱內界面波傳播的作用力機制的研究仍不夠完善。有鑒于此，本文在Wei等[32]研究的基礎上研究重力突降后貯箱體系內各相似準則數之間的關聯。

1 計算模型

圖1顯示了空間推進劑貯箱的常用構型，其關于z軸呈旋轉對稱，中間圓柱段高為0.8 m，半徑為1.5 m；上下兩端封頭為橢球形，其短軸長度為1.5 m。貯箱上部增壓氣體為氦氣，下部推進劑為液氧，充液率為50%。表1給出了相應的流體物性參數，對應溫度78 K、壓力2×105Pa、重力加速度g方向沿z軸負方向。

表1 氣液兩相流體物性參數Table 1 Material parameters of the gas and liquid phases

圖1 空間推進劑貯箱構型Fig.1 Space propellant tank configuration

采用流體體積法（VOF）[33]求解液氧貯箱內氣液自由界面在變重力環境下的運動。數值計算中假定流動為不可壓縮黏性層流流動，相應的控制方程如下：

其中，

其中，V為流體運動速度，p為流體壓力，Fσ為采用連續表面張力模型[34]計算表面張力作用所對應的等效體積力。ρm為總的流體密度，μm為總的流體動力學黏性系數，α為液體的體積分數函數，1-α為氣體的體積分數函數，下標2表示氣相。

貯箱壁面處流體運動滿足無滑移條件，即有如下表示：

其中，Vw為貯箱壁面的運動速度，這里取為0。

壓力p和體積分數函數α滿足Neumann邊界條件，沿壁面法向梯度為零，即

接觸角恒定，在三相接觸線位置表示如下：

其中，n為液氣相界面處的單位法向量，nw和tw分別為壁面的單位法向量和單位切向量，θ為固定接觸角。

初始時貯箱內流體靜止，初始壓力設為2×105Pa，且自由液面保持水平，這意味著初始狀態對應于毛細效應可忽略的強重力情況。數值模擬初始時刻，重力階躍下降至設定的殘余重力水平。本文數值模擬中，Bo數范圍為1～5000，其大小是通過改變殘余重力水平來實現的，特征長度取為貯箱柱段半徑R。

本研究中的物理對象本身具有對稱性，故選用軸對稱二維模型進行數值計算。計算域網格通過ICEM來生成，近壁區適當細化（見圖1）。控制方程中的非定常項采用一階隱式格式，對流項采用QUICK格式，擴散項采用Least Squares Cell Based格式，壓力離散采用PRESTO算法，VOF采用Geo-Reconstruct格式離散。壓力-速度耦合方程的求解采用SIMPLE算法。

為了提升求解精度，同時也為了能夠較清晰地捕捉到更大Bo數的較為細小的波動，計算中選擇網格數量為36萬。計算模型的有效性驗證及網格無關性檢驗可參照Wei等[32]的工作。

2 結果分析與討論

圖2給出了重力突降（Bo=200）后不同時刻界面位形波的演化過程，更詳細的描述參見文獻[32]的描述。圖2中用虛線箭頭標注了近壁波谷（本文將其稱為初始波谷）的運動，其運動特征可以表征重力突降引起的初始界面波運動規律。

圖2 Bo=200時初始界面波的移動（虛線箭頭表示近壁波谷）Fig.2 Displacement of the main wave trough at Bo = 200 (Dashed arrow indicates near-wall troughs)

重力突降之后，毛細作用使得近壁區液體爬升，進而導致鄰近區域液面下沉，對自由界面形成擾動。這一擾動將導致界面波的形成，并由近壁區域向貯箱中心傳播。明確可見界面波的出現需要一定時間，可將近壁下沉區域后端波峰（即高過初始液面的隆起）的出現時刻作為界面波形成的初始時刻ti。另外，氣液界面波總是以波群的形式向外傳播，波群中速度最快的波會率先在貯箱中心匯合產生波反射，影響后續界面波的傳播。將首個波傳播到貯箱軸線位置的時刻定義為界面波自由傳播的終止時刻te。將上述兩個時刻之間的過程作為界面波自由傳播階段。

圖3給出了重力突降后不同Bo數條件下貯箱內氣液界面上初始波谷徑向位置r隨時間的變化。可以看到界面波自由傳播階段初始波谷的位移呈現出明顯的線性特征，其斜率即界面波傳播速度隨Bo數增大而加大。

圖3 不同Bo數下初始波谷位置的變化Fig.3 Movement of the initial trough of interface wave under different values of the Bond number

當Bo數大于1時，重力突降之后，貯箱內界面波傳播的特征速度可表示為[35]

設定物性恒定，則Bo數越大，重力突降之后的殘余重力水平越高，由式（9）可知，此時貯箱內的界面波的特征傳播速度自然隨之增大。為了進一步考察界面波傳播速度的具體變化規律，定義無量綱速度U′為界面波真實傳播速度u和特征速度u0之比，即

圖4顯示了無量綱速度U′隨Bo數的變化情況。在所研究的Bo數范圍內，無量綱波速的變化趨勢呈現出明顯的非線性變化特征。隨著Bo數的增大，無量綱波速趨向于一個定值1.414，這對應著極高重力條件下重力主導界面波傳播的極限情況。

圖4 不同Bo數下界面波的無量綱傳播速度變化情況Fig.4 Change of the dimensionless propagation velocity of the interface wave under different values of the Bond number

重力突降后，毛細作用將克服重力引起液體沿壁面的爬升。毛細爬升的響應時間可表示如下[35]：

據此可定義界面波自由傳播階段的起始與終止的無量綱時間分別為

圖5給出了不同Bo數下界面波自由傳播的無量綱終止時刻。可以看到，當Bo≤38時，界面波自由傳播階段的無量綱終止時刻小于1，此時界面波的運動完全由壁面附近的毛細力爬升所驅動，阻礙液面爬升的是液體慣性，液面毛細力和液體慣性力間的平衡關系即

圖5 不同Bo數下界面波自由傳播的無量綱終止時刻Fig.5 Dimensionless termination moments of free propagation of interface waves for different Bo numbers

從式（14）可以得到，此時界面波傳播的特征速度U～(σ/ρR)1/2。引入Froude數Fr=u/（gR）1/2表示慣性力與重力的相對比值，可以得到Bo數較小時界面波傳播速度滿足如下標度關系：

擬合本文數值模擬結果可以得到比例系數為3.143，即Fr=3.143Bo-1/2。利用式（9）可以將Fr數改寫成如下形式：

聯立式（15）和式（16），可以得到

考慮到在小擾動波傳播過程中，小波長 （λ?lc，其中lc= （σ/ρg）1/2為毛細長度或Laplace長度）情形滿足U′～（λ/lc）-1/2。鑒于λ/lc= （λ/R）Bo1/2，式（17）表明在小Bo數條件下，無量綱擾動波長λ/R將是一個不依賴于Bo數的常數。基于數值模擬結果，可以得到λ/R≈ 0.1。

圖6給出了不同Bo數下界面波自由傳播的無量綱初始時刻。可以看到，當Bo≥1800時，界面波自由傳播階段的無量綱初始時刻將大于1。這表明近壁毛細主導的爬升已經停止，液體慣性作用消退，界面波動特征（例如波長）將由重力和毛細力決定，即毛細重力波。相應地，界面波的波數k可以表示為如下形式[36]：

圖6 不同Bo數下界面波自由傳播的無量綱初始時刻Fig.6 Dimensionless initial moments of free propagation of interface waves for different Bo numbers

毛細重力波的群速度可表示為

結合式（18），可得

根據圖4中大Bo數極限U′趨近于1.414，上式中“～”可直接代之以“=”，進一步可以得到如下標度關系：

對于中等Bo數，界面波的運動由壁面處的毛細爬升和毛細重力波的傳播二者共同主導，即毛細力、慣性力和重力共同主導界面波的傳播，相關過程特征可以采用Churchill-Usagi關聯方法進行擬合。為了更好表征過渡區特性，引入衰減效應，對經典的Churchill-Usagi關聯方法進行了修正，得到如下標度關系：

式中，A，C和n為待定擬合參數，依賴于過渡區實際特征。基于本文數值模擬結果的最小二乘法擬合，得到如下可涵蓋Bo數全區域的標度關系：

圖7給出了式（23）與數值模擬結果的比較，可以看出二者吻合得很好。此外，基于數值模擬結果與大、小Bo數極限情形中變化趨勢的偏離程度，建議在重力突降情形中，氣液界面波自由傳播過程中的毛細爬升區終止于Bo=20，而毛細重力波區則起始于Bo=2000。

圖7 Fr數與Bo數之間的標度關系Fig.7 Scaling relationship between the Fr number and Bo numbers

3 結論

數值模擬了貯箱內因重力突降引起的氣液界面波的自由傳播過程，研究了Bo數從1～5000范圍內界面波的傳播規律，得到如下的主要結論。

（1）貯箱內界面波的傳播速度隨著Bo數的增加而增大。

（2）當Bo≥2000時，界面波的傳播由重力和毛細力共同決定，處于毛細重力波區，其特征傳播速度為(gσ/ρ)1/4，傳播速度滿足標度規律Fr=1.414Bo-1/2。

（3）當Bo≤20時，界面波的傳播由毛細力和慣性力共同作用的毛細爬升所決定，處于毛細爬升區，其特征傳播速度為(σ/ρR)1/2，傳播速度滿足標度規律Fr=3.143Bo-1/4。

（4）采用引入衰減效應的Churchill-Usagi關聯方法進行擬合，得到了可以涵蓋Bo數全范圍的界面波傳播速度標度關系，統一描述了毛細爬升區、過渡區和毛細重力波區界面波自由傳播過程特征。