微重力下圓管毛細(xì)流動解析近似解研究*

李永強(qiáng) 張晨輝 劉玲 段俐 康琦

1)(東北大學(xué)理學(xué)院應(yīng)用力學(xué)研究所,沈陽 110819)

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

(2012年9月3日收到;2012年9月8日收到修改稿)

1 引言

目前關(guān)于微重力環(huán)境下的圓管毛細(xì)流動問題受到了廣泛研究,并得到了不同的液面高度對時(shí)間的依賴關(guān)系.Lucas[1]和Washburn[2]考慮了液面的毛細(xì)力和液體在水平管中的摩擦力,利用力平衡方程得到了毛細(xì)管中液面高度h隨時(shí)間的變化關(guān)系為h-t1/2,并且Lucas-Washburn方程被Bell和Cameron[3]的實(shí)驗(yàn)所證實(shí).由于Lucas-Washburn模型在液體流動的初始階段假定液面速度為無窮大,因此它的有效性受到了限制,并且實(shí)驗(yàn)[2?5]都是在地面重力條件下用很小尺寸的管徑完成的,因而在實(shí)驗(yàn)中觀察到了液體的Lucas-Washburn行為(h-t1/2).Siegel[6]首次利用自由落體平臺進(jìn)行了微重力條件下的實(shí)驗(yàn),微重力時(shí)間為0.74 s,在實(shí)驗(yàn)中觀測到了液面高度隨時(shí)間變化的線性階段(h-t),但是在進(jìn)入Lucas-Washburn階段之前實(shí)驗(yàn)就結(jié)束了.Petrash等[7]首次利用落塔實(shí)驗(yàn)研究了液面高度隨時(shí)間的變化,該實(shí)驗(yàn)中微重力時(shí)間為2.25 s,其結(jié)果也表明存在液面高度與時(shí)間呈線性關(guān)系的階段,除此以外,他們還發(fā)現(xiàn),在進(jìn)入線性階段之前,還應(yīng)有一個(gè)初始的h-t2關(guān)系.有研究者嘗試通過導(dǎo)入幾種其他的影響,像慣性力[8,9]、液池自由面的曲率[7]、動態(tài)接觸角[6,10,11]等來改進(jìn)Lucas-Washburn方程.Levine等[12]提出了目前為止關(guān)于毛細(xì)流動的最詳盡的理論.Stange等[13]利用德國不萊梅的落塔研究了微重力條件下圓柱形管中液體的毛細(xì)流動過程,從理論上建立了圓管毛細(xì)流動的非線性動力學(xué)方程,通過數(shù)值解和實(shí)驗(yàn)相結(jié)合的方法給出了毛細(xì)流動過程在不同階段時(shí)毛細(xì)爬升高度隨時(shí)間變化的規(guī)律.Wang等[14]研究了微重力環(huán)境下插入液體中的圓管毛細(xì)流動過程,研究了接觸角和管徑尺寸對毛細(xì)流動的影響,實(shí)驗(yàn)結(jié)果表明:隨著接觸角的增加,圓管中的液面高度h和液面速度˙h隨著微重力時(shí)間單調(diào)的減小.隨著管徑尺寸的增加,在一定的微重力時(shí)間內(nèi),h-t曲線及˙h-t曲線并不是單調(diào)變化的,不同管徑尺寸的h-t曲線可能相互發(fā)生交叉.

目前關(guān)于圓管毛細(xì)流動非線性動力學(xué)方程的求解尚無詳盡的解析近似解研究,均采用數(shù)值解的方法.關(guān)于非線性動力學(xué)方程的解析近似求解常采用攝動法,如龐加萊-林茲斯泰德法、漸進(jìn)法、多尺度法和平均法等,但這些方法對弱非線性系統(tǒng)是有效的,而對強(qiáng)非線性系統(tǒng)則難以應(yīng)用.廖世俊[15]基于代數(shù)拓?fù)鋵W(xué)中連續(xù)映射的概念,提出了一種新的求解非線性問題解析解的方法—–同倫分析法.不同于所有已知的解析近似方法,同倫分析法通過引入輔助參數(shù)和輔助函數(shù)來調(diào)節(jié)和控制級數(shù)解的收斂區(qū)域和收斂速度,該方法對強(qiáng)、弱非線性系統(tǒng)都適合,已成功解決工程技術(shù)中的許多非線性問題[16?25].

本文采用同倫分析法求解了Stange等建立的圓管毛細(xì)流動非線性動力學(xué)方程,以級數(shù)解的形式獲得了內(nèi)角流動的近似解析解,并給出了級數(shù)解的具體表達(dá)形式.

2 同倫分析法基本思想

Stange等[13]從理論上對微重力條件下插入液面的圓管毛細(xì)流動進(jìn)行了分析,幾何構(gòu)型的示意圖如圖1所示[26].在容器內(nèi)壁與圓管外壁設(shè)有防爬擋板,圓管插入液面的深度記為l0,圓管內(nèi)液面中心點(diǎn)的位置用h表示.

圖1 圓管毛細(xì)流動幾何構(gòu)型的示意圖

初始液體高度記為l0(插入液面深度),圓管內(nèi)液面高度為h,液面速度為˙h.圓管的半徑為R,液面曲率的主半徑為R1和R2,動態(tài)接觸角為γd.容器內(nèi)液體自由面的半徑Rc可由防爬擋板間的距離a以及容器中自由面的中心線半徑c計(jì)算得出.Stange等[13]在Levine等[12]方法的基礎(chǔ)上做了改進(jìn),提出了毛細(xì)流動的非線性動力學(xué)方程為

其中

式中Re是雷諾數(shù),Rex=x/ν,Red=d/ν分別是基于流動長度和流動直徑d的雷諾數(shù).一般情況下,可取K(x)≈4/3.方程(1)的邊界條件為

Stange等[13]對三種介質(zhì)和七種容器的圓管毛細(xì)流動進(jìn)行了實(shí)驗(yàn)和數(shù)值分析,介質(zhì)的特性參數(shù)如表1所示,容器結(jié)構(gòu)參數(shù)如表2所示.其中表1中的SF 0.65和SF 1.00分別指運(yùn)動黏度為0.65和1.00 cSt(1 cSt=10?4m2/s)的硅油,FC-77是含碳的氟化液系列中的一種.

表1 實(shí)驗(yàn)介質(zhì)參數(shù)

表2 實(shí)驗(yàn)容器結(jié)構(gòu)參數(shù)

由表 1和表 2可知,當(dāng)采用 SF 1.00為實(shí)驗(yàn)介質(zhì)和第Ⅰ號容器時(shí),tr達(dá)到最小值,此時(shí)trmin=0.008127 s;當(dāng)采用FC-77為實(shí)驗(yàn)介質(zhì)和第Ⅶ號容器時(shí),tr達(dá)到最大值,此時(shí)trmax=0.94024 s.當(dāng)tr分別達(dá)到最大值和最小值時(shí),s(t)的變化曲線如圖2所示.

圖2 s(t)-t曲線

由圖2可知,s([t)≈1,故本文中設(shè)]s(t)=1.將方程(1)中的進(jìn)行泰勒展開,經(jīng)計(jì)算發(fā)現(xiàn)取前四階即可滿足計(jì)算精度,則

其中b1,b2,b3,b4和b5為展開系數(shù),則方程(1)可簡化為

其中

按照廖世俊[15]同倫分析法的基本思想,根據(jù)方程(5)的形式,定義非線性算子

其中,q∈[0,1]為嵌入變量,Φ(t;q)是未知函數(shù)h(t)的映射.設(shè)零階形變方程為

該方程滿足邊界條件

當(dāng)q從0變化到1時(shí),Φ(t;q)從初始值h0(t)變化到(5)式的解h(t).

根據(jù)(5)式,選取線性輔助算子

該線性算子具有性質(zhì)

其中,C為任意常數(shù).根據(jù)邊界條件(3)式,初始值h0(t)設(shè)為

其中α為未知常數(shù).利用泰勒展開定理,將Φ(t;q)展開成關(guān)于q的冪級數(shù)

其中

如果輔助參數(shù)ˉh、輔助函數(shù)H(t)和未知常數(shù)α選取合適的值,則當(dāng)q=1時(shí),(13)式能夠收斂,則可得

(15)式必定是非線性方程(5)的某一個(gè)解.定義矢量

對前面求得的零階形變方程(7)求輔助參數(shù)q的n階導(dǎo)數(shù),再令q=0,然后除以n!,可得n階形變方程

該方程滿足邊界條件

其中

根據(jù)解表達(dá)原則和系數(shù)遍歷性原則[15]可以令輔助函數(shù)H(t)為

其中k為整數(shù),研究中發(fā)現(xiàn)k≤2時(shí),hn(t)中含有l(wèi)n(t)項(xiàng),為了避免出現(xiàn)ln(t)項(xiàng),取

通過求解方程(17),發(fā)現(xiàn)hn(t)可以表達(dá)成

其中bn,j為系數(shù),將(23)式帶入到方程(16)中,按照t的次冪相等原則就可得到系數(shù)bn,j的具體表達(dá)形式.例如當(dāng)n=1時(shí),可得

系數(shù)bn,j推導(dǎo)的相關(guān)工作可由符號運(yùn)算軟件Maple或Mathmatica進(jìn)行.則當(dāng)級數(shù)取前m階時(shí),可得

式中

其中INT[]表示取整函數(shù).

3 計(jì)算結(jié)果分析

通過上述分析可知,僅需選擇合適的輔助參數(shù)α和ˉh時(shí),可確保級數(shù)解(15)式收斂.將級數(shù)解(15)式帶入到方程(1)中,取平方后在區(qū)域t∈[0,5]上積分,可得平方殘差為

其中

圖3為初始液體高度l0=50 mm,三種實(shí)驗(yàn)介質(zhì)在編號為Ⅶ(即R=35.0 mm)的容器內(nèi)毛細(xì)流動時(shí),液面高度取前4階級數(shù)解時(shí),方程(1)的平方殘差E和α,ˉh的關(guān)系曲線圖.

圖3 三種實(shí)驗(yàn)介質(zhì)在容器Ⅶ內(nèi),液面高度前4階級數(shù)解的平方殘差Em與ˉh和α的關(guān)系曲線(a)SF 0.65;(b)SF 1.00;(c)FC-77

表3 不同級數(shù)下,(24)式中系數(shù)ck的值(實(shí)驗(yàn)介質(zhì):FC-77,容器編號:Ⅶ,l0=50 mm)

從圖3(a)可以看出,Em的最小值在α∈[0.02,0.05],ˉh∈[?0.005,0.005]范圍內(nèi),級數(shù)解在該區(qū)域內(nèi)收斂;從圖3(b)可以看出,Em的最小值在α∈[?2.5×10?4,2.5×10?4],ˉh∈[0.02,0.06]范圍內(nèi);從圖3(c)可以看出,Em的最小值在α∈[?0.001,0.001],ˉh∈[?0.01,?0.005]范圍內(nèi).為了精確獲得最小平方殘差位置,方程(27)的極值點(diǎn)為

聯(lián)立求解方程(29)即可得到α和ˉh的值.

當(dāng)實(shí)驗(yàn)介質(zhì)為FC-77,容器編號為Ⅶ(即R=35.0 mm),初始液體高度l0=50 mm時(shí),不同級數(shù)下(26)式中系數(shù)ck的值如表3所示.

實(shí)驗(yàn)介質(zhì)為FC-77,容器編號為Ⅶ時(shí),液面高度h(t)在不同時(shí)刻下的m階同倫近似解與龍格庫塔(R-K)數(shù)值計(jì)算結(jié)果的對比如表4所示.從表4可以看出,當(dāng)級數(shù)m=4時(shí),級數(shù)解就開始收斂,所以在圓管毛細(xì)流動動力學(xué)方程計(jì)算中,取級數(shù)m=4.

表4 h(t)的m階同倫近似解與R-K數(shù)值計(jì)算結(jié)果在不同時(shí)刻的對比(實(shí)驗(yàn)介質(zhì):FC-77,容器編號:Ⅶ,l0=50 mm)

圖4為取級數(shù)m=4時(shí),初始液體高度l0=50 mm,三種實(shí)驗(yàn)介質(zhì)在編號為Ⅰ(即R=2.0 mm)的容器內(nèi)毛細(xì)流動時(shí),液面高度h(t)的同倫解析近似解(HAM)與R-K數(shù)值解對比圖.由圖4可知,應(yīng)用同倫分析法得到的解析近似解與數(shù)值法求得的解是相當(dāng)符合的.

圖4 h(t)的前4階同倫解析近似解曲線

實(shí)驗(yàn)介質(zhì)為FC-77時(shí),(26)式中的系數(shù)ck隨初始液體高度l0的變化情況如圖5所示,由于系數(shù)較多,這里僅給出了c1的變化情況,其他兩種實(shí)驗(yàn)介質(zhì)的c1變化規(guī)律與FC-77相同.從圖5可以看出,c1隨l0和圓管半徑R的增加迅速減小,當(dāng)圓管半徑R＞10 mm時(shí),c1接近于零.如果實(shí)驗(yàn)時(shí)想要獲得較高的液面高度h(t),則圓管半徑不宜大于10 mm.

圖5 實(shí)驗(yàn)介質(zhì)FC-77在不同容器時(shí),系數(shù)c1與初始液體高度l0的關(guān)系曲線

4 結(jié)論

本文應(yīng)用同倫分析方法獲得了微重力環(huán)境下圓管毛細(xì)流動的解析近似解并給出了級數(shù)解的表達(dá)形式.同倫分析方法為我們提供了一個(gè)方便的方法來控制漸進(jìn)級數(shù)的收斂,這是同倫分析方法和其他方法根本性的區(qū)別.從同倫分析法與四階R-K法的計(jì)算結(jié)果比較表明,同倫分析法具有較好的計(jì)算精度.

[1]Lucas V R 1918 Kolloid-Z.23 15

[2]Washburn E W 1921 Phys.Rev.17 273

[3]Bell JM,Cameron F K 1906 J.Phys.Chem.10 658

[4]Rideal E K 1922 Philos.Mag.44 1152

[5]LeGrand E J,Rense WA 1945 J.Appl.Phys.16 843

[6]Siegel R 1961 J.Appl.Mech.83 165

[7]Petrash D A,Nelson T M,Otto E W 1963 NASA TN D-1582

[8]Jeje A A 1979 J.Colloid Interf.Sci.69 420

[9]Ichikawa N,Satoda Y 1994 J.Colloid Interf.Sci.162 350

[10]Joos P,Remoortere P,Bracke M 1990 J.Colloid Interf.Sci.136 189

[11]Qu′er′e D 1997 Europhys.Lett.39 533

[12]Levine S,Reed P,Watson E J,Neale G 1976 In Colloid and Interface Science(New York:Academic)p403

[13]Stange M,Dreyer M E,Rath H J 2003 Phys.Fluids 15 2587

[14]Wang C X,Xu S H,Sun Z W,Hu WR 2009 AIAA J.11 2642

[15]Liao S J 2006 Beyond Perturbation:Introduction to the Homotopy Analysis Method(Beijing:Science Press)p204(in Chinese)[廖世俊2006超越攝動—–同倫分析方法導(dǎo)論(北京:科學(xué)出版社)第204頁]

[16]Cheng J,Liao S J 2007 Acta Mech.Sin.39 715(in Chinese)[成均,廖世俊2007力學(xué)學(xué)報(bào)39 715]

[17]Liao S J 2003 J.Fluid Mech.488 189

[18]Li Y Q,Zhu D W,Li F 2009 Chin.J.Mech.Eng.45 37(in Chinese)[李永強(qiáng),朱大巍,李鋒2009機(jī)械工程學(xué)報(bào)45 37]

[19]Li Y Q,Li F,Zhu D W 2010 Compos.Struct.92 1110

[20]Yuan P X,Li Y Q 2010 Appl.Math.Mech.31 1293

[21]Li Y Q,Li L,He Y L 2011 Compos.Struct.93 360

[22]Li Y Q,Zhu D W 2011 Compos.Struct.93 880

[23]Shi Y R,Yang H J 2010 Acta Phys.Sin.59 67(in Chinese)[石玉仁,楊紅娟2010物理學(xué)報(bào)59 67]

[24]Yang P,Chen Y,Li Z B 2010 Acta Phys.Sin.59 3668(in Chinese)[楊沛,陳勇,李志斌2010物理學(xué)報(bào)59 3668]

[25]Liao S J 2012 Homotopy Analysis Method for Nonlinear Differential Equations(Beijing:Higher Education Press)p285

[26]Dreyer M E 2007 Spring Tracts in Mordern Physics 221 51

[27]Sparrow E M,Lin S H,Lundgren T S 1964 Phys.Fluids 7 338