微重力條件下部分充液貯箱氣液界面波動特性的數(shù)值模擬1)

魏列 杜王芳 趙建福 李凱

(中國科學(xué)院力學(xué)研究所微重力重點(diǎn)實(shí)驗(yàn)室,北京 100190)

(中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

引言

空間航行出現(xiàn)以來,航天器的復(fù)雜性以及任務(wù)的艱巨性不斷提高,推進(jìn)劑管理系統(tǒng)對空間環(huán)境的適應(yīng)性要求也越來越嚴(yán)苛[1-2].推進(jìn)劑貯箱中液體行為的天地差異主要來自由于體積力的減弱或消失,控制流體行為的毛細(xì)力作用發(fā)生了顯著改變[3].液體火箭在飛行過程中必須考慮氣液界面如何響應(yīng)重力水平改變的問題,同時確保在微重力滑行之后具備可靠的引擎再啟動能力.微重力環(huán)境下,火箭上面級貯箱中氣液界面的運(yùn)動及氣液兩相流動與分布特性一直都是重要的研究內(nèi)容,也是先進(jìn)空間流體管理技術(shù)研發(fā)的關(guān)鍵[4].

為了能夠在地面常重力環(huán)境研究微重力條件下的流體界面現(xiàn)象,Masica 和Petrasb[5]用劉易斯研究中心的落塔研究了圓柱容器里氣液界面的重定位運(yùn)動,發(fā)現(xiàn)重定位過程中界面中心點(diǎn)向液相運(yùn)動的速率是恒定的,其大小可由經(jīng)驗(yàn)關(guān)系式給出.Labus 和Masica[6]同樣利用劉易斯研究中心的落塔給出了球形容器中自由界面的重定位過程.依據(jù)實(shí)驗(yàn)結(jié)果,他們推測液面涌泉的產(chǎn)生取決于收集界面處液體的流速.Dreyer等[7-9]利用不萊梅落塔研究重力水平變小后自由界面的運(yùn)動及界面振蕩的阻尼特性.Jenson等[10]利用國際空間站實(shí)驗(yàn)研究了固定接觸線和自由接觸線邊界條件對容器內(nèi)自由液面受到加速度擾動時的運(yùn)動特征的影響.魏延明和潘海林[11]用落塔實(shí)驗(yàn)探索了表面張力貯箱中推進(jìn)劑的靜平衡和重定位過程中的液面變化.張景芳[12]用落塔實(shí)驗(yàn)研究了球錐貯箱內(nèi)自由液面在較寬邦德數(shù)范圍內(nèi)的準(zhǔn)靜態(tài)平衡位形.顧方[13]利用落塔實(shí)驗(yàn)研究了微重力條件下有無側(cè)向干擾對半球形貯箱內(nèi)液體重定位過程的影響.莊保堂等[14]用落塔實(shí)驗(yàn)研究了板式表面張力貯箱內(nèi)蓄液器的蓄留特性.Liu等[15-16]利用落塔實(shí)驗(yàn)研究了不同充液比和貯箱方位對葉片式貯箱內(nèi)推進(jìn)劑的界面位形及爬升速度的影響.Li等[17]用落塔實(shí)驗(yàn)和數(shù)值模擬研究了微重力條件下水作為工質(zhì)的自由液面的振蕩頻率與貯箱尺寸和充液比的關(guān)系.

如上所述,地基落塔短時微重力實(shí)驗(yàn)一直以來在空間流體管理的研究和設(shè)計(jì)中占據(jù)著重要地位,但微重力持續(xù)時間短是其最大缺陷.數(shù)值模擬能夠自主加載邊界條件,且數(shù)值模擬所計(jì)算的時間范圍更加自由.李章國等[18-19]用流體體積法來捕捉貯箱內(nèi)液體推進(jìn)劑在不同重力加速度和接觸角下氣液界面的形位分布及變化.姜志杰等[20]采用“Surface Evolver”軟件研究了充液比、接觸角、邦德數(shù)對球形貯箱內(nèi)自由液面分布的影響.趙婷婷等[21]對微重力條件下推進(jìn)劑管理裝置模型的界面位形進(jìn)行了數(shù)值模擬,分析了界面特征點(diǎn)的速度與壓力的變化.Li等[22-23]采用動態(tài)接觸角模型數(shù)值研究了貯箱內(nèi)自由界面在不同充液比以及不同重力水平下的振動特性及自由界面的阻尼振蕩特性.胡齊等[24]數(shù)值仿真了板式表面張力貯箱內(nèi)推進(jìn)劑在不同重力加速度下的重定位過程.黃濱等[25]用數(shù)值模擬方法研究了側(cè)向加速度環(huán)境下板式表面張力貯箱內(nèi)推進(jìn)劑在不同充液比時的自由液面演化規(guī)律.嚴(yán)天等[26]用數(shù)值模擬方法對微重力條件下板式貯箱內(nèi)低溫推進(jìn)劑在初始液位傾斜、側(cè)向和旋轉(zhuǎn)激勵作用下流體的運(yùn)動特性開展了瞬態(tài)仿真.陳磊等[27]采用流體體積法兩相流動模型對微重力條件下板式貯箱內(nèi)流體行為以及加注過程中推進(jìn)劑的流動特性進(jìn)行數(shù)值仿真.周振華等[28]用數(shù)值仿真方法研究了微重力條件下儲液器內(nèi)氣液界面狀態(tài)與儲液器尺寸、壁面浸潤性和工質(zhì)充液比之間的關(guān)聯(lián)性.肖立明等[4]數(shù)值研究了貯箱內(nèi)自由液面的演化特征,但未能對界面波動現(xiàn)象的相似特性進(jìn)行分析.

以上研究對理解貯箱內(nèi)的流體運(yùn)動特性有諸多幫助,但關(guān)于界面波傳播規(guī)律及相關(guān)流動現(xiàn)象的相似準(zhǔn)則等的定量表征尚不明確.有鑒于此,本文基于微重力流體物理最常采用的邦德數(shù)相似準(zhǔn)則,針對空間貯箱常用構(gòu)型和典型尺寸,設(shè)計(jì)了3 個縮比模型,采用流體體積法方法和連續(xù)表面張力模型,數(shù)值模擬了原型貯箱和縮比模型中因加速度變化引起的貯箱內(nèi)氣液兩相流動及氣液界面上界面波的傳播過程,以期加深對界面波的傳播規(guī)律和流動相似性的理解與掌握.

1 計(jì)算模型

1.1 物理模型

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

圖1 原型貯箱示意圖Fig.1 Schematic diagram of the prototype tank

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

1.2 數(shù)學(xué)模型和計(jì)算方法

本文采用流體體積法[29]求解液氧貯箱內(nèi)氣液自由界面在變重力環(huán)境下的運(yùn)動.

本文計(jì)算中假定流動為不可壓縮黏性層流流動,相應(yīng)的控制方程分別表示如下

式(1)～ 式(3)中,V為流體運(yùn)動速度,ρ為流體密度,p是流體壓力,μ是流體動力學(xué)黏性系數(shù),Fσ是采用連續(xù)表面張力模型計(jì)算的氣液界面附近表面張力作用對應(yīng)的等效體積力,f是流體體積函數(shù),為任意空間單元內(nèi)流體所占體積與該單元可容納流體體積之比.

貯箱壁面處流體運(yùn)動滿足無滑移條件,即

其中,Vw是貯箱壁面的運(yùn)動速度,這里取為0.

流體體積函數(shù)在壁面附近滿足接觸角條件,本文不考慮接觸角滯后效應(yīng),假設(shè)表觀接觸角恒定,即

初始時流體靜止,液面保持水平,這意味著初始狀態(tài)對應(yīng)于毛細(xì)效應(yīng)可忽略的強(qiáng)重力情況.

利用ANSYS Fluent CFD 來求解控制方程.本研究中的物理對象本身具有對稱性,故選用軸對稱二維模型進(jìn)行數(shù)值計(jì)算.計(jì)算域網(wǎng)格通過ICEM 來生成,近壁區(qū)適當(dāng)細(xì)化(見圖1).控制方程中的非定常項(xiàng)采用一階隱式格式,對流項(xiàng)采用QUICK 格式,擴(kuò)散項(xiàng)采用least squares cell based 格式,壓力離散采用PRESTO 算法,VOF 采用Geo-Reconstruct 格式離散.壓力-速度耦合方程的求解采用SIMPLE 算法.

2 基于邦德數(shù)相似準(zhǔn)則的縮比模型

縮比模型試驗(yàn)是在地面上研究微重力條件下流體行為的重要試驗(yàn)方法[30].在研究液體火箭和宇宙飛船的過程中,有很多與推進(jìn)劑運(yùn)動特性有關(guān)的問題需要用模擬試驗(yàn)的方法來解決[31].根據(jù)相似理論,模型試驗(yàn)的相似準(zhǔn)則數(shù)要和原型的同名相似準(zhǔn)則數(shù)相等.一般來說,不可能讓原型和模型之間所有的相似準(zhǔn)則數(shù)都相同,對于模型試驗(yàn),應(yīng)盡可能使得原型和模型之間主要的相似準(zhǔn)則數(shù)相等.

鑒于地面重力環(huán)境浮力的主導(dǎo)作用和空間微重力環(huán)境浮力作用被抑制后表面張力作用將凸顯出來,微重力流體物理中邦德數(shù)(Bo=Δρa(bǔ)L2/σ) 自然成為一個主要的相似準(zhǔn)則數(shù),表征了重力和表面張力的相對強(qiáng)弱.考慮原型貯箱在空間微重力環(huán)境(約為10-5g0,這里g0表示地面重力水平) 邦德數(shù)Bo=140,根據(jù)邦德數(shù)相似準(zhǔn)則,本文設(shè)計(jì)了3 組縮比模型(表2),來研究原型與縮比模型中的流動相似特征.

表2 縮比模型參數(shù)設(shè)計(jì)Table 2 Parameter design of the scale-down models

在重力突變后,接觸線區(qū)域的運(yùn)動導(dǎo)致了氣液界面產(chǎn)生波動,并沿界面?zhèn)鞑?鑒于Bo=140 ＞ 1,本文研究中的界面波屬于毛細(xì)重力波,其特征速度[3]可表示為

這樣,界面波傳播到容器中心的特征時間是

利用式(6)和式(7),可以得到如下標(biāo)度關(guān)系

式中,t和u分別為重力突變后貯箱內(nèi)液體所經(jīng)歷的時間和速度大小.

3 網(wǎng)格無關(guān)性和模擬有效性驗(yàn)證

3.1 網(wǎng)格無關(guān)性檢驗(yàn)

本文通過CFD-POST 軟件處理每一個計(jì)算時間步的數(shù)值仿真結(jié)果,利用其內(nèi)置的功能去定位幅度最大的波谷及其鄰的波峰,給出波谷隨時間變化的坐標(biāo)值以求得界面波的傳播速度.由于并不存在有限強(qiáng)度的單色波,波速較大的波會先到達(dá)中心匯合產(chǎn)生波的反射,干擾主波的行進(jìn).本文選取中心點(diǎn)第一個波谷形成的時刻作為干擾的起始.此外,受限于網(wǎng)格分辨率和近壁擾動,主波波峰形成的初始時間需要依據(jù)波傳播特征確定.由這兩個時間節(jié)點(diǎn)所確定時段屬于界面波自由傳播階段,為壁面兩側(cè)擾動對數(shù)據(jù)的影響,取其中部約80%部分通過對波谷位置變化的分析確定界面波的傳播速度.

圖2 顯示了原型貯箱和模型Ⅲ的表面張力波波速隨網(wǎng)格數(shù)的變化.可以看到,隨著網(wǎng)格加密、網(wǎng)格數(shù)增大,界面波速度變化逐漸減小,計(jì)算結(jié)果趨于收斂;當(dāng)網(wǎng)格數(shù)增大到一定程度,繼續(xù)增加網(wǎng)格數(shù)對計(jì)算結(jié)果改進(jìn)不大.綜合考慮計(jì)算效率和精度,在本文后續(xù)分析中,原型貯箱網(wǎng)格數(shù)選為16 萬,模型Ⅲ網(wǎng)格數(shù)選為24 萬.模型Ⅰ,Ⅱ網(wǎng)格數(shù)分別為16 萬和24 萬.

圖2 網(wǎng)格無關(guān)性檢驗(yàn)Fig.2 Mesh independency

3.2 模擬有效性驗(yàn)證

Concus[32]給出了直立圓柱中靜液面的漸近解析解,其中,界面爬升高度的解析表達(dá)式如下

其中

式中,B表示邦德數(shù),θ 表示接觸角.液體的爬升高度h隨不同容器圓柱段半徑的變化及其數(shù)值模擬與漸近解的對比如表3 所示.爬升高度是通過壁面接觸點(diǎn)和中間平直段的高度差的平均求得.可以看出,爬升高度的數(shù)值解和漸近解吻合得很好.

由于液體黏度很小,數(shù)值計(jì)算需要相當(dāng)長的時間才能到達(dá)一個相對穩(wěn)定的漸近穩(wěn)態(tài)界面,因此,對數(shù)值平衡界面的提取比較困難,提取的漸近穩(wěn)態(tài)界面并不是完全光滑的.本文數(shù)值模擬的漸近穩(wěn)態(tài)界面和Concus[32]漸近解析解的對比如圖3 所示,其中的無量綱長度都是以相應(yīng)的圓柱段半徑為特征長度得到的.可以看出,數(shù)值模擬計(jì)算出來的平衡位形與理論上的液面靜平衡是完全吻合的.這與表3 爬升高度的對比在有效性驗(yàn)證上是相互補(bǔ)充的.

圖3 平衡液面位形數(shù)值解和解析解的對比Fig.3 Comparison of numerical solution and asymptotic solution of equilibrium liquid surface configuration

表3 液體爬升高度的漸近解和數(shù)值解對比Table 3 Asymptotic and numerical simulation results of liquid climbing height

圖4 中把原型貯箱內(nèi)自由面運(yùn)動的仿真結(jié)果與肖立明等[4]仿真結(jié)果進(jìn)行了比較,可以看到兩者的界面演化結(jié)果在相同時刻下是基本一致的.此外,本文計(jì)算的原型貯箱中表面張力波傳播速度為1.19 cm/s,這與肖立明等[4]相似條件下的數(shù)值模擬結(jié)果1.25 cm/s 基本一致,二者相對差值約5%.

圖4 原型貯箱內(nèi)自由面的運(yùn)動Fig.4 Numerical simulations of free surface motion in the prototype tank

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

4.1 氣液界面形態(tài)

圖5 中給出了原型和縮比模型在同一無量綱時間下界面波的運(yùn)動過程,這里的無量綱時間定義與式(7)相同.

圖5 原型和模型自由面運(yùn)動仿真對比Fig.5 Comparison of free surface motion simulations between prototype and models

可以看到,在氣泡形成之前,原型和模型的氣液界面的運(yùn)動基本上是一樣的.壁面處的三相接觸線在無量綱時間為0.28 時爬升到最大高度,隨后緩慢下移.在無量綱時間為0.80 的時候,主波在容器中心融合,形成一個小氣泡,在中心處還可以看到出現(xiàn)了一個很小的涌泉.

不過,界面波在貯箱中軸線匯聚形成的氣泡形狀略有不同.原型中氣泡形狀近似為三角形,且體積最大;而在3 個縮比模型中,氣泡較小且形狀較為光滑.其原因可能和網(wǎng)格尺寸所對應(yīng)的物理空間尺度等因素有關(guān):由于縮比模型的物理尺寸減小,在計(jì)算網(wǎng)格數(shù)不變甚至增多的情況下,縮比模型的網(wǎng)格分辨率總是更高的.盡管略有差異,縮比模型與原型貯箱內(nèi)部的氣液兩相界面變化過程之間的相似性還是極為明顯和確定的.

4.2 界面波速與流場速度

表4 列出了原型和模型計(jì)算結(jié)果,其中,界面波速uw為上述界面波自由傳播時段的傳播速度,其數(shù)值在界面波自由傳播階段是恒定的;氣泡初始形成時間t0即圖5 中的無量綱時間t*=0.80,對應(yīng)于界面波首次在中軸線處的匯聚和中心處氣泡的形成.表中還列出了韋伯?dāng)?shù)We=ρRuw2/σ及其相對于原型數(shù)值的相對誤差.

表4 表明界面波速隨著縮比增大(即物理尺寸減小)而增大,但對應(yīng)的韋伯?dāng)?shù)近似恒定,最大相對誤差不到2%.根據(jù)弗勞德數(shù)定義Fr=uw2/(gR),邦德數(shù)恒定與韋伯?dāng)?shù)近似恒定可以直接推論出弗勞德數(shù)同樣近似恒定.這說明在目前的邦德數(shù)相似準(zhǔn)則約束下,原型貯箱和縮比模型之間近似存在韋伯?dāng)?shù)相似準(zhǔn)則,或等價地,弗勞德數(shù)相似準(zhǔn)則.這說明在相關(guān)體系中,重力、表面張力和慣性力的作用,均保持了相對恒定的強(qiáng)弱程度.

表4 波速分析Table 4 Wave speed analysis

圖6 是不同模型中某個對應(yīng)位置處韋伯?dāng)?shù)隨無量綱時間的演化過程.顯然,韋伯?dāng)?shù)相似準(zhǔn)則不僅對貯箱內(nèi)界面波傳播過程適用,對任意位置的氣液兩相流動也是適用的.這說明界面波速可以用來表征貯箱內(nèi)流體的流動特征.

圖6 不同模型流場對應(yīng)點(diǎn)處Weber 數(shù)的演化特征(圖中選擇的對應(yīng)點(diǎn)的無量綱坐標(biāo)為(0.60,-0.46))Fig.6 The evolution of the Weber number at the corresponding locations in different models (the dimensionless coordinates of the corresponding points selected in the figure are (0.60,-0.46))

圖6 也顯示出不同縮比模型中流動特性在若干時刻存在細(xì)微的差異,特別是局部流動的韋伯?dāng)?shù)或局部速度的極大值會隨著模型縮比的增大而有所提高.表4 界面波速也顯示出同樣的變化趨勢.

產(chǎn)生上述差異的原因,一方面源于縮比模型物理尺寸減小導(dǎo)致的表面張力作用(以氣液界面兩側(cè)毛細(xì)壓力差2σ/R來表征)增強(qiáng),使得重力突變后所能釋放出來的驅(qū)動力增強(qiáng),流速峰值增加.另一方面,體系物性不變時,根據(jù)近似的韋伯?dāng)?shù)相似準(zhǔn)則,隨著模型縮比比例增大(即貯箱尺寸減小),系統(tǒng)內(nèi)的流動速度將增大.大的速度和小的物理尺寸必然導(dǎo)致貯箱內(nèi)流體剪切率增大,強(qiáng)化了黏性耗散作用,這與苗楠等[30]的結(jié)果是一致的.

因此,在空間大尺度貯箱的地面常重力模擬試驗(yàn)設(shè)計(jì)中,需要對縮比比例進(jìn)行適當(dāng)限制,以免產(chǎn)生過大的縮比誤差.

5 結(jié)論

本文采用數(shù)值模擬方法對微重力環(huán)境下部分充液貯箱內(nèi)氣液界面和氣、液兩相介質(zhì)的運(yùn)動特征進(jìn)行了研究.針對空間貯箱常用構(gòu)型和典型尺寸,基于微重力流體物理領(lǐng)域常用的邦德數(shù)相似準(zhǔn)則設(shè)計(jì)了3 個縮比模型,比較了原型貯箱和縮比模型中因重力突變引起的貯箱內(nèi)氣液兩相流動特性及氣液界面上界面波的傳播特征.對于大型空間推進(jìn)劑貯箱,即使在微重力環(huán)境,殘余重力依然使得系統(tǒng)邦德數(shù)明顯大于1,貯箱內(nèi)氣液兩相體系將同時受重力、表面張力以及黏性力等作用.數(shù)值模擬結(jié)果驗(yàn)證了基于邦德數(shù)相似準(zhǔn)則設(shè)計(jì)的縮比模型與原型貯箱之間的運(yùn)動相似性,發(fā)現(xiàn)在滿足邦德數(shù)相似準(zhǔn)則的前提下,系統(tǒng)還近似滿足韋伯?dāng)?shù)相似準(zhǔn)則,或等價地,近似滿足弗勞德數(shù)相似準(zhǔn)則.此外,數(shù)值模擬結(jié)果也表明原型貯箱和縮比模型間的運(yùn)動存在細(xì)微偏差.根據(jù)近似的韋伯?dāng)?shù)相似準(zhǔn)則可知,縮比比例加大,貯箱尺寸減小,重力突變后由表面張力釋放出來的驅(qū)動力增強(qiáng),相同韋伯?dāng)?shù)下流動速度將增大,黏性耗散作用隨之增強(qiáng).

本文研究結(jié)果可以用于指導(dǎo)空間貯箱流體管理技術(shù)地面模擬試驗(yàn)的方案設(shè)計(jì)等.