李永強 張晨輝 劉玲 段俐 康琦

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

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

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

1 引言

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

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

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

2 同倫分析法基本思想

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

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

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

其中

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

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

表1 實驗介質(zhì)參數(shù)

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

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

圖2 s(t)-t曲線

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

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

其中

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

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

該方程滿足邊界條件

當(dāng)q從0變化到1時,Φ(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時,(13)式能夠收斂,則可得

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

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

該方程滿足邊界條件

其中

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

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

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

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

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

式中

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

3 計算結(jié)果分析

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

其中

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

圖3 三種實驗介質(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的值(實驗介質(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)的極值點為

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

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

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

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

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

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

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

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

4 結(jié)論

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

[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é)報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)[李永強,朱大巍,李鋒2009機械工程學(xué)報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é)報59 67]

[24]Yang P,Chen Y,Li Z B 2010 Acta Phys.Sin.59 3668(in Chinese)[楊沛,陳勇,李志斌2010物理學(xué)報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