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

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

1)(東北大學理學院應用力學研究所,沈陽 110819)

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

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

1 引言

目前關于微重力環境下的圓管毛細流動問題受到了廣泛研究,并得到了不同的液面高度對時間的依賴關系.Lucas[1]和Washburn[2]考慮了液面的毛細力和液體在水平管中的摩擦力,利用力平衡方程得到了毛細管中液面高度h隨時間的變化關系為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階段之前實驗就結束了.Petrash等[7]首次利用落塔實驗研究了液面高度隨時間的變化,該實驗中微重力時間為2.25 s,其結果也表明存在液面高度與時間呈線性關系的階段,除此以外,他們還發現,在進入線性階段之前,還應有一個初始的h-t2關系.有研究者嘗試通過導入幾種其他的影響,像慣性力[8,9]、液池自由面的曲率[7]、動態接觸角[6,10,11]等來改進Lucas-Washburn方程.Levine等[12]提出了目前為止關于毛細流動的最詳盡的理論.Stange等[13]利用德國不萊梅的落塔研究了微重力條件下圓柱形管中液體的毛細流動過程,從理論上建立了圓管毛細流動的非線性動力學方程,通過數值解和實驗相結合的方法給出了毛細流動過程在不同階段時毛細爬升高度隨時間變化的規律.Wang等[14]研究了微重力環境下插入液體中的圓管毛細流動過程,研究了接觸角和管徑尺寸對毛細流動的影響,實驗結果表明:隨著接觸角的增加,圓管中的液面高度h和液面速度˙h隨著微重力時間單調的減小.隨著管徑尺寸的增加,在一定的微重力時間內,h-t曲線及˙h-t曲線并不是單調變化的,不同管徑尺寸的h-t曲線可能相互發生交叉.

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

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

2 同倫分析法基本思想

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

圖1 圓管毛細流動幾何構型的示意圖

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

其中

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

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

表1 實驗介質參數

表2 實驗容器結構參數

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

圖2 s(t)-t曲線

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

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

其中

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

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

該方程滿足邊界條件

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

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

該線性算子具有性質

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

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

其中

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

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

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

該方程滿足邊界條件

其中

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

其中k為整數,研究中發現k≤2時,hn(t)中含有ln(t)項,為了避免出現ln(t)項,取

通過求解方程(17),發現hn(t)可以表達成

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

系數bn,j推導的相關工作可由符號運算軟件Maple或Mathmatica進行.則當級數取前m階時,可得

式中

其中INT[]表示取整函數.

3 計算結果分析

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

其中

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

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

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

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

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

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

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

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

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

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

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

圖5 實驗介質FC-77在不同容器時,系數c1與初始液體高度l0的關系曲線

4 結論

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

[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超越攝動—–同倫分析方法導論(北京:科學出版社)第204頁]

[16]Cheng J,Liao S J 2007 Acta Mech.Sin.39 715(in Chinese)[成均,廖世俊2007力學學報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機械工程學報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物理學報59 67]

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