陳上通 吳 笛 王 佳 段 俐 康 琦,

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

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

引言

表面張力驅(qū)動(dòng)流,在此定義為由表面張力引起的附加壓力驅(qū)動(dòng)的自發(fā)界面流,是航天器貯箱內(nèi)液體行為的重要組成部分.為了在微重力環(huán)境下下高效地管理液體,有必要深入研究表面張力問(wèn)題.

Washburn[1]首次提出了毛細(xì)管內(nèi)液體爬升的運(yùn)動(dòng)方程.他們證實(shí)該過(guò)程可以用毛細(xì)驅(qū)動(dòng)壓力、黏滯阻力和重力之間的平衡來(lái)描述,推導(dǎo)出了著名的Lucus-Washburn 方程.Concus 和Finn[2]探索了具有內(nèi)角的容器中的液體平衡界面,并提出了著名的Concus-Finn 條件.Levine 等[3]考慮了表面張力驅(qū)動(dòng)壓力、黏滯阻力和對(duì)流損失,推導(dǎo)出了圓管內(nèi)液體爬升高度的二階微分方程.Stange 等[4]提出了一個(gè)更全面的圓管內(nèi)液體爬升模型,它補(bǔ)充考慮了彎液面的重定位過(guò)程、動(dòng)態(tài)接觸角和泊肅葉流動(dòng)的啟動(dòng)過(guò)程.Jiang 等[5]在1979 年改進(jìn)了動(dòng)態(tài)接觸的計(jì)算方法并被廣泛接受.Dreyer 等[6]探究了平行平板之間的表面張力驅(qū)動(dòng)流動(dòng)并得到了液體爬升高度方程.Weislogel 和Lichter[7]獲得了液體沿尖內(nèi)角處的流動(dòng)方程,并擴(kuò)展到由不同潤(rùn)濕性形成的尖內(nèi)角處[8]和圓內(nèi)角處[9].Yue 和Wang[10]采用數(shù)值方法研究了三維自由表面的動(dòng)力學(xué)問(wèn)題.Higuera 等[11]研究了潤(rùn)濕液體在兩個(gè)垂直板之間形成一個(gè)小角度的狹窄間隙中的流動(dòng)過(guò)程.Wolf 等[12]提出了一種基Lattice-Boltzmann 方法、通過(guò)考慮流體和固壁之間長(zhǎng)程相互作用的影響來(lái)模擬平行平板之間的液體爬升.Bolleddula 等[13]提出了一種新的沿不同截面形狀尖內(nèi)角處流動(dòng)的解析解,并探討了實(shí)驗(yàn)?zāi)Ｐ蛶缀涡螤畹挠绊?魏月興等[14]研究了內(nèi)角處表面張力驅(qū)動(dòng)流動(dòng)理論并將其應(yīng)用到板式貯箱的結(jié)構(gòu)設(shè)計(jì)上.李永等[15-17]利用落塔試驗(yàn)和數(shù)值模擬探究了板式貯箱的流體管理性能,得到了豐富的微重力試驗(yàn)數(shù)據(jù).Reyssat[18]探究了平板和圓弧板之間形成的液橋界面形貌特征,建立的模型可以準(zhǔn)確預(yù)測(cè)液橋外形.Wu 等[19]研究了彎曲尖內(nèi)角處的流動(dòng)過(guò)程,通過(guò)修正尖內(nèi)角曲率半徑帶來(lái)的影響,得到了流動(dòng)距離與曲率半徑的關(guān)系.Romain 等[20]研究了軸對(duì)稱幾何模型中的表面張力驅(qū)動(dòng)流動(dòng).Chen 等[21]利用毛細(xì)驅(qū)動(dòng)流動(dòng)理論對(duì)貯箱內(nèi)PMD 結(jié)構(gòu)進(jìn)行了優(yōu)化,并開(kāi)展了落塔實(shí)驗(yàn)和基于VOF 法的仿真分析,檢驗(yàn)了其性能.Cheng 等[22]研究了變直徑圓管中的流動(dòng)問(wèn)題.Chen等[23]研究了微重力下橢圓管內(nèi)的表面張力驅(qū)動(dòng)流動(dòng),并提出了兩段式的流動(dòng)新模型.

本文探究了微重力環(huán)境下成一定角度平板間的表面張力驅(qū)動(dòng)流,考慮了動(dòng)態(tài)接觸角、對(duì)流壓力損失、黏滯阻力和液池內(nèi)彎曲自由面的影響,并利用基于VOF 方法的數(shù)值模擬進(jìn)行驗(yàn)證.此外,微分方程可變換成由作用在板間控制體上一系列合力的方程,通過(guò)同時(shí)考慮兩個(gè)主要作用力,將整個(gè)爬升過(guò)程分成三個(gè)階段.

1 理論推導(dǎo)

研究模型如圖1 所示,包含圓柱形液池和自上方插入液池一定深度的平板模型,液池內(nèi)壁和平板外壁均有防爬擋板,選擇直角坐標(biāo)系來(lái)分析該問(wèn)題.初始液面為水平面,z 向原點(diǎn)位于水平面上.液體爬升高度為h,液體爬升平均速度為.平板水平截面如圖2 所示.平板左端的距離是2a,右端的距離是2c,板的寬度是b.平板上特定位置與和x 軸間的距離是d.液池內(nèi)彎曲界面的曲率半徑是用防爬擋板之間的距離i 與液池內(nèi)氣液界面中心線和平板模型對(duì)稱線之間的距離j 計(jì)算.為了計(jì)算i 與j 的值,平板之間的梯形入口等效成半徑為re的圓形入口,半徑re為理論分析的基本假設(shè)是

圖1 研究模型正向剖視圖Fig.1 Front view of the model

圖2 平板模型的水平截面Fig.2 Horizontal cross section of the plates

(1)氣液界面上的剪切力忽略不計(jì);

(2)流動(dòng)過(guò)程是恒溫的;

(3)流動(dòng)是充分發(fā)展的泊肅葉流;

(4)流體是牛頓流體,不可壓縮且均質(zhì);

(5)液體與壁面之間無(wú)滑移.

平板上特定位置與和x 軸間的距離為

流動(dòng)是充分發(fā)展的泊肅葉流,液體僅具有z 向速度分量u,流動(dòng)邊界條件為

u 為關(guān)于x,y,t 的函數(shù),t 為時(shí)間.可以認(rèn)為板間區(qū)域由x 向上的無(wú)數(shù)個(gè)dx 分段組成,在每個(gè)dx 分段內(nèi)u 為關(guān)于y,t 的函數(shù),與x 無(wú)關(guān).結(jié)合式(2)并根據(jù)N-S 方程可得每個(gè)dx 分段內(nèi)速度場(chǎng)分布為

結(jié)合平板間流量公式

式中 Ω 代表板間水平截面全域.可將式(3)轉(zhuǎn)化為

N-S 方程z 向分量去除零項(xiàng)后為

式中 ρ 為密度,ν 為液體運(yùn)動(dòng)學(xué)黏性系數(shù).對(duì)該式在板間水平截面上對(duì)x 和y 進(jìn)行積分,得到

將式(7)對(duì)z 從z=0 到z=h 進(jìn)行積分,可以得到

采用與文獻(xiàn)[14]相似的處理方法,可以得到板間頂端凹液面的壓力為

式中 p0為液面附近大氣壓,σ 為表面張力,αd和 αs分別為動(dòng)態(tài)接觸角和靜態(tài)接觸角.對(duì)于完全潤(rùn)濕的液體,即靜態(tài)接觸角為0°,其動(dòng)態(tài)接觸角有一經(jīng)驗(yàn)公式

式中μ 為動(dòng)力學(xué)黏性系數(shù).該計(jì)算模型受到了廣泛運(yùn)用和驗(yàn)證,文獻(xiàn)[4,6,20]均采用該計(jì)算模型.

為了計(jì)算板間入口處的作用力,在等效圓形入口附近建立半徑為re的半球形控制體2,如圖3所示.

圖3 板間入口處的等效半球形控制體Fig.3 Equivalent circular entrance and the control volume around the inlet

式中 pR為液池內(nèi)彎曲界面引起的附加壓力,為

液池中自由面的曲率半徑Rc的近似表達(dá)式為

在板間入口處的作用力沿z 向分力為

經(jīng)半球面R=re進(jìn)入控制體2 的流量z 向分量為

經(jīng)板間入口流出控制體2 的流量z 向分量為

經(jīng)半球面R=re進(jìn)入控制體2 的加速度流z 向分量為

經(jīng)入口處流出控制體2 的加速度流z 向分量為

因此,半球形控制體2 內(nèi)的整體加速度的近似表達(dá)式為

他的腳離我的頭只有幾英寸遠(yuǎn)，我理應(yīng)去安慰他，我本應(yīng)該主動(dòng)去安慰他才對(duì)，因?yàn)槲覐男【褪鞘苓@種教育長(zhǎng)大的。相反，我覺(jué)得那樣做很惡心。看起來(lái)那么強(qiáng)壯的人，不應(yīng)該表現(xiàn)得這么脆弱。為什么他不能像其余人一樣悄悄地哭呢？

結(jié)合式(11)～ 式(18)及控制體2 上的動(dòng)量方程,可以求得板間入口處的作用力z 向分量為

將式(1)、式(5)、式(9)、式(10)和式(20)代入式(8),可以得到液體爬升高度h 的二階微分方程

該二階微分方程可利用四階Runge-Kutta 法并結(jié)合初始條件 h(t=0)=h˙(t=0)=0 求解.

2 數(shù)值模擬

本文建立的三維計(jì)算模型如圖4(a)所示.由于液池內(nèi)氣液界面彎曲造成的影響很小,為了建模的方便,選擇邊長(zhǎng)為140 mm 的方形液池來(lái)代替理論分析中的圓柱形液池,并在理論計(jì)算中使用79 mm的等效半徑,該簡(jiǎn)化所造成的影響可以忽略不計(jì),模型的總高度為120 mm.建立的網(wǎng)格模型如圖4(b)所示,網(wǎng)格總數(shù)約為100 萬(wàn),在所有壁面附近建立邊界層網(wǎng)格,首層邊界層的高度約為0.1 mm,相鄰兩層網(wǎng)格膨脹比為1.2.由于計(jì)算模型簡(jiǎn)單,網(wǎng)格質(zhì)量好,調(diào)整邊界層網(wǎng)格后進(jìn)行第二次仿真,兩次結(jié)果差異很小,驗(yàn)證了網(wǎng)格對(duì)仿真結(jié)果的無(wú)關(guān)性.同時(shí)為了減小仿真可能產(chǎn)生的隨機(jī)誤差,取兩次結(jié)果的平均值作為最終結(jié)果.

圖4 仿真模型Fig.4 Numerical model

本文使用KF-96 系列硅油開(kāi)展研究工作,該系列型號(hào)硅油以其運(yùn)動(dòng)學(xué)黏性系數(shù)作為標(biāo)號(hào),比如2 號(hào)硅油(SF 2)、5 號(hào)硅油(SF 5)和10 號(hào)硅油(SF 10),他們的物性參數(shù)如表1 所示.本文選擇層流作為流動(dòng)模式,壓力-速度耦合方程由SIMPLIC 算法進(jìn)行數(shù)值求解,壓力空間離散式方程選用“Body Force Weighted”算法求解,梯度空間離散化方程選擇“Least Square Cell”算法求解,動(dòng)量空間離散化方程使用二階迎風(fēng)格式求解,松弛因子采用默認(rèn)設(shè)置,當(dāng)方程迭代余量減少到10-6時(shí),計(jì)算被認(rèn)為收斂.在計(jì)算過(guò)程中,Courant 數(shù)(Co)大部分時(shí)候小于1,這表明計(jì)算過(guò)程非常穩(wěn)定.Co 對(duì)于瞬態(tài)流動(dòng)很重要.對(duì)于一維網(wǎng)格,它的定義為

表1 硅油物性參數(shù)(25 °C)Table 1 Liquid properties (25 °C)

式中,u 代表流速,Δx 代表網(wǎng)格尺寸,Δ t 代表時(shí)間步長(zhǎng).

板間流動(dòng)的仿真?zhèn)纫晥D如圖5 所示,黃色曲面代表氣液界面.在t=0 s 時(shí),液體全部位于液池底部,氣液界面為一平面.仿真開(kāi)始后,液體快速向上爬升并在板間形成一彎曲液面.液體與平板的接觸線不是直線,如圖6 所示.因此,為了計(jì)算液體爬升高度,選擇接觸線的兩端點(diǎn)H1,H3和最高點(diǎn)H2,將其平均值(H1+H2+H3)/3 視為液體爬升高度h.

圖5 仿真結(jié)果側(cè)視圖.液體為10 號(hào)硅油,a=2 mm,b=14 mm,c=5 mmFig.5 Front view of the numerical model as the liquid is SF 10 and a=2 mm,b=14 mm,c=5 mm

圖6 板上液相分布,紅色部分代表液體Fig.6 Liquid distribution in a plate and the red part represents liquid

理論和仿真爬升高度的對(duì)比如圖7 所示.本文使用了6 種尺寸的平板模型.圖中實(shí)線代表理論結(jié)果,方塊代表仿真結(jié)果,并用不同顏色代表不同型號(hào)硅油的爬升數(shù)據(jù).黑色代表2 號(hào)硅油的數(shù)據(jù),紅色代表5 號(hào)硅油的數(shù)據(jù),藍(lán)色代表10 號(hào)硅油的數(shù)據(jù).由于模型高度的限制,對(duì)于2 號(hào)硅油,取其2 s 內(nèi)的爬升數(shù)據(jù);對(duì)于5 號(hào)硅油,取其4 s 內(nèi)的爬升數(shù)據(jù);對(duì)于10 號(hào)硅油,取其6 s 內(nèi)的爬升數(shù)據(jù).本文忽略了板間液面從一開(kāi)始的平面到彎曲平衡界面的發(fā)展過(guò)程.此外,對(duì)于爬升高度的測(cè)量方法也將帶來(lái)一定誤差.總體來(lái)看,數(shù)值結(jié)果與理論結(jié)果吻合良好.對(duì)比圖8(a)和圖8(b),或者圖8(c)和圖8(d)可以看出,相同黏性的硅油,在較寬間距板間的爬升速度比較小間距板間的慢,平板之間的距離越寬,板間流速越慢.從理論分析中可知,表面張力的增加有助于增加毛細(xì)驅(qū)動(dòng)壓力,加快液體流速,而液體黏性的增加會(huì)減小流速.

圖7 液體爬升高度隨時(shí)間變化Fig.7 Liquid climbing height h vs.time t

3 機(jī)理分析

為了進(jìn)一步分析,可以將式(14)轉(zhuǎn)化成作用在控制體1 上合力的方程,如下

式中各力的含義為

(1)板間毛細(xì)驅(qū)動(dòng)壓力

(2)液池內(nèi)慣性力

(3)板間入口對(duì)流壓力損失

(4)板間黏滯阻力

(5)板間慣性力

(6)液池內(nèi)黏滯阻力

(7)液池內(nèi)附加壓力

不同工況下各力發(fā)展情況如圖8 所示,橫坐標(biāo)是時(shí)間t,縱坐標(biāo)是力F.深紅色、綠色、淡藍(lán)色、深藍(lán)色和洋紅色實(shí)線分別代表平板間毛細(xì)驅(qū)動(dòng)壓力Fcp、液池中的慣性力Fir、板間入口處的對(duì)流壓力損失Fpl、兩塊平板上的黏滯阻力Ffp和液池內(nèi)的附加壓力Fcr,黑色虛線代表平板間的慣性力Fip.毛細(xì)驅(qū)動(dòng)壓力從最大值開(kāi)始變化,迅速減小至最小值后又開(kāi)始緩慢增加,該變化是由于動(dòng)態(tài)接觸角的變化.從上述四圖中可以看出,無(wú)論是高Oh數(shù)(Oh=還是低Oh數(shù)工況,液池內(nèi)慣性力、板間入口壓力損失和板上黏滯阻力均先后發(fā)揮了主要作用,這點(diǎn)與橢圓管和同心圓管內(nèi)的流動(dòng)特征大不相同.液池內(nèi)慣性力從最大值開(kāi)始變化,迅速減小,直至被對(duì)流壓力損失超越;入口處對(duì)流壓力損失從0 開(kāi)始變化,迅速增加到最大值后開(kāi)始緩慢減小.最終在某一時(shí)刻,板上黏滯阻力超越入口壓力損失發(fā)揮主要作用.因此,板間的表面張力驅(qū)動(dòng)流動(dòng)均可以劃分成三個(gè)階段.在第一階段中,令毛細(xì)驅(qū)動(dòng)壓力等于液池內(nèi)慣性力,可以得到

圖8 各力隨時(shí)間變化情況Fig.8 Forces’ development vs time

該方程的初始條件為 h˙(t=0)=0 .結(jié)合初始條件可求得第一階段的爬升高度為

為了計(jì)算第二階段的速度方程,令毛細(xì)驅(qū)動(dòng)壓力等于入口處壓力損失,可以得到

令第一、第二階段的速度方程相等,可以求出劃分這兩個(gè)階段的時(shí)刻t1

結(jié)合第二階段的初始條件h1(t1),可得到第二階段爬升高度為

為了計(jì)算第三階段的速度方程,令毛細(xì)驅(qū)動(dòng)壓力等于板上黏滯阻力,可以得到

令第二、第三階段的速度方程相等,可以求出劃分這兩個(gè)階段的時(shí)刻t2為

結(jié)合第三階段的初始條件h2(t2),可得第三階段的爬升高度為

不同工況的計(jì)算參數(shù)如表2 所示,利用Oh 數(shù)區(qū)分不同工況.

相比于高Oh 數(shù)工況,低Oh 數(shù)工況下第二階段持續(xù)時(shí)間更長(zhǎng).對(duì)比圖8(a)和圖8(b),或者圖8(c)和圖8(d),均可看出在低Oh 數(shù)時(shí)入口對(duì)流壓力損失發(fā)揮主要作用的時(shí)間更長(zhǎng),這意味著第二階段持續(xù)時(shí)間更長(zhǎng).例如,在第四個(gè)仿真工況中,其流動(dòng)Oh 數(shù)為0.004 50,第二階段持續(xù)0.647 s;而在第六個(gè)仿真工況中,其流動(dòng)Oh 數(shù)為0.022 1,第二階段僅持續(xù)0.039 s,遠(yuǎn)小于低Oh 數(shù)工況下的持續(xù)時(shí)間.楔形平板間的表面張力驅(qū)動(dòng)流可以分成三個(gè)階段,如圖9(a)所示,x 軸是時(shí)間t,y 軸是液體爬升高度h,藍(lán)線代表爬升高度與時(shí)間的關(guān)系,圖中標(biāo)注了三個(gè)階段的范圍.第一階段中爬升曲線為拋物線型,液體爬升高度與t2成正比;第二階段中爬升曲線為直線型,液體爬升高度與t 成正比;第三階段液體爬升高度與t1/2成正比.理論和仿真爬升高度的對(duì)比如圖9(b)所示,其中黑色實(shí)現(xiàn)代表理論結(jié)果,方塊代表仿真結(jié)果,劃分不同流動(dòng)階段的時(shí)刻t1和t2均標(biāo)注在圖中,為了更直觀地觀察爬升高度的變化規(guī)律,建立了一條從時(shí)刻t1開(kāi)始的虛直線.仿真結(jié)果和理論結(jié)果吻合良好,爬升高度的變化規(guī)律也與理論分析一致,充分驗(yàn)證了理論分析的正確性.

4 結(jié)論

本文通過(guò)理論分析,推導(dǎo)出了微重力環(huán)境下成一定夾角平板間表面張力驅(qū)動(dòng)流爬升高度的二階微分方程,并利用基于VOF 方法的數(shù)值模擬開(kāi)展了驗(yàn)證工作.本文采用了6 種研究模型和3 種型號(hào)硅油,理論結(jié)果與仿真結(jié)果吻合良好.

本文通過(guò)進(jìn)一步的理論分析,還發(fā)現(xiàn)同時(shí)考慮兩個(gè)主要作用力,可將平板間的表面張力驅(qū)動(dòng)流動(dòng)分為三個(gè)階段,在這三個(gè)階段中液體爬升高度h 分別與t2,t 和t1/2成正比.在第一階段,液池中的慣性力和毛細(xì)驅(qū)動(dòng)壓力發(fā)揮主要作用.但是液池中的慣性力迅速減小,流動(dòng)進(jìn)入第二階段,由對(duì)流壓力損失和毛細(xì)驅(qū)動(dòng)壓力發(fā)揮主要作用.并且,在高Oh 數(shù)情況下,第二階段的持續(xù)時(shí)間很短暫.第三階段時(shí)由板上黏滯阻力和毛細(xì)驅(qū)動(dòng)壓力發(fā)揮主要作用.相較于第三階段,第一、第二階段持續(xù)時(shí)間都很短暫,整個(gè)流動(dòng)過(guò)程以第三階段為主.本章研究結(jié)果可為板式貯箱中導(dǎo)流板和蓄液葉片等的設(shè)計(jì)和空間流體管理提供理論依據(jù).