重力跳變引起的貯箱內(nèi)氣液界面波傳播規(guī)律*

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

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

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

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

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

0 引 言

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

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

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

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

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

1 計算模型

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

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

圖1 空間推進劑貯箱構(gòu)型Fig.1 Space propellant tank configuration

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

其中，

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

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

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

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

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

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

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

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

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

2 結(jié)果分析與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

聯(lián)立式（15）和式（16），可以得到

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

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

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

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

結(jié)合式（18），可得

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

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

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

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

圖7 Fr數(shù)與Bo數(shù)之間的標度關(guān)系Fig.7 Scaling relationship between the Fr number and Bo numbers

3 結(jié)論

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

（1）貯箱內(nèi)界面波的傳播速度隨著Bo數(shù)的增加而增大。

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

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

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