不同重力環(huán)境和傾斜角度對毛細管內(nèi)界面振蕩行為的影響*

霍曉智 王 晴 顧君蘋 王戰(zhàn)濤 于 強 王慶功

1(中國科學(xué)院國家空間科學(xué)中心 北京 100190)

2(中國科學(xué)院大學(xué)計算機科學(xué)與技術(shù)學(xué)院 北京 100049)

3(中國空間技術(shù)研究院 錢學(xué)森空間技術(shù)實驗室 北京 100094)

4(中山職業(yè)技術(shù)學(xué)院 機電工程學(xué)院 中山 528400)

0 引言

毛細現(xiàn)象對人類生活及空間活動均有重要的意義。在微重力條件下，由于重力的缺失，毛細力占主導(dǎo)地位，由毛細力驅(qū)動的流體分布是空間熱控、流體管理的關(guān)鍵基礎(chǔ)[1-3]。毛細結(jié)構(gòu)內(nèi)的流體動力學(xué)行為受毛細尺度、流體物性、壁面潤濕性、毛細管傾斜角度等因素的影響[4-7]，在重力變化時會誘導(dǎo)界面發(fā)生振蕩現(xiàn)象，毛細結(jié)構(gòu)內(nèi)兩相界面的運動過程較為復(fù)雜。

20世紀初，Lucas-Washburn方程（L-W方程）的提出從理論上建立了重力條件下毛細管內(nèi)界面上升的規(guī)律，并證明了毛細界面上升前期高度h與時間t的關(guān)系（h2～t）。然而，L-W方程忽略了慣性對毛細動力學(xué)的影響。針對這一問題，Bosanquet[8]將液柱的動量變化加入L-W方程中，Szekely等[9]、Levine等[10]從進入毛細管口流體的能量入手修正方程，從而更加準確地描述了毛細管內(nèi)界面運動演變規(guī)律。Quere[11]通過實驗觀察到毛細管內(nèi)界面的振蕩過程，Wang等[5]通過理論分析研究了毛細管內(nèi)界面振蕩機理。

隨著研究的不斷深入，對常規(guī)毛細現(xiàn)象已經(jīng)有較清楚的認識。然而，在空間微重力環(huán)境下，毛細管內(nèi)的界面爬升過程有所變化。隨著重力的消失，界面振蕩行為也隨之消失。Dreyer等[12]、Stange等[13]在落塔中進行了微重力環(huán)境不同尺度的毛細實驗，提出毛細界面運動可分為三個階段：初始階段以毛細力為驅(qū)動力，第二階段以入口處流體流動引起的對流力為主導(dǎo)力，在第三階段，管間摩擦力對毛細運動的影響較大。Wang 等[2]研究了微重力環(huán)境毛細管尺度對毛細驅(qū)動行為的影響，發(fā)現(xiàn)隨著毛細管尺寸的增加，毛細界面運動進入最后一個階段的時間也隨之增加。此外，Xu等[14]、Chen等[15,16]研究了不同形狀毛細管對毛細管內(nèi)界面爬升動力學(xué)的影響。

在空間及地面應(yīng)用中，毛細管并非一直處于與重力加速度垂直的方向。由于毛細管的傾斜或重力環(huán)境的變化，毛細管內(nèi)的界面過程會受到進一步影響。Gao等[7]通過實驗研究了傾斜角度對毛細管液柱上升高度的影響，所得實驗結(jié)果與L-W方程存在明顯差異。Tian等[17]和Barroso等[18]通過理論建模的方式研究了Taylor-Hauksbee Cell傾斜時的毛細上升問題，證明了高度與時間遵循h(huán)～t1/3的規(guī)律。

目前，對于非常重力環(huán)境和傾斜毛細管內(nèi)的界面振蕩現(xiàn)象理論分析仍不充分。為解析這一問題，本文建立了傾斜毛細管內(nèi)界面變化過程的理論模型，分析了毛細界面振蕩的動力學(xué)效應(yīng)。同時針對毛細管建立數(shù)值計算模型，開展了模擬研究，獲得了可視化結(jié)果。將理論分析、數(shù)值模擬結(jié)果和現(xiàn)有實驗結(jié)果進行對比，證明了兩類模型的可靠性，同時詳細分析了不同重力環(huán)境和傾斜角度對毛細管界面上升及振蕩行為的影響規(guī)律。

1 理論模型

1.1 傾斜毛細管內(nèi)界面位置隨時間的演化規(guī)律

如圖1所示，毛細現(xiàn)象是液體依靠毛細力作用在毛細管內(nèi)自動爬升的過程。依據(jù)能量平衡原理，毛細上升和/或下降的動態(tài)過程不僅與表面能）、動能）、重力勢能）、凈輸入能）、黏性能的耗散（）有關(guān)，還與入口/出口處壓力損失（） 的瞬態(tài)變化有關(guān)[5]。因此，能量守恒定律遵循以下公式：

圖1 傾斜毛細管內(nèi)彎液面位置與結(jié)構(gòu)Fig.1 Schematic of meniscus when capillary tube is tilted

表面能的變化率為

其中，R，Δz，l的 表示如圖1所示，σ為表面張力。

動能的變化率為

其中，ρ表示液體密度，v為速度。

重力勢能的變化率為

其中，g為重力加速度，α為毛細管相對于垂直方向的傾斜角度。

質(zhì)量變化率為˙，表達式如下：

流入毛細管中的凈動能變化率為

流入毛細管中的凈重力勢能變化率為

由式(6)和式(7)可得毛細管中凈輸入能變化率為

根據(jù)Hagen-Poiseuille關(guān)系，黏性阻力產(chǎn)生的壓降可表示為

其中，μ為液體的黏度，A為截面積；毛細管中產(chǎn)生的黏性阻力為

由黏性阻力產(chǎn)生的黏性耗散能量變化率可寫為

根據(jù)流體力學(xué)理論，壓頭損失可以計算為

其中，ξ為液體運動的局部阻力系數(shù)，可通過液池（A1）和毛細管（A2）的橫截面積之比估算（見圖2）[5]，即

圖2 壓頭損失計算Fig.2 Schematic diagram of head pressure loss

結(jié)合上述公式，液體上升/下降的能量平衡方程可寫為

當(dāng)毛細管傾斜角度α=0時，毛細管垂直于液面，式(15)與文獻[5]一致。在一般情況中，液體對固體的潤濕情況為部分潤濕狀態(tài)。因此，在計算該類問題時，需要考慮接觸角的變化，即在式(15)的左邊乘以接觸角的余弦值（cosθ）。

1.2 動態(tài)接觸角

液體與毛細管壁的接觸角變化主要是由三相線處受力不平衡引起的，即

其中，θd表示液體與毛細管壁的動態(tài)接觸角。

毛細力（Fc）為

其中，θ0為液體靜止時與毛細管壁的接觸角。

摩擦力（Ff）為

其中，?f為接觸線摩擦因數(shù)，cpin為釘扎力，可以用接觸角滯后量化，u為接觸線處的速度。

結(jié)合式(16)～(18)可得毛細過程中的動態(tài)接觸角，即

則液體界面上升/下降的能量平衡方程如下：

1.3 Runge-Kutta求解微分方程

為求解上述微分方程，采用四階Runge-Kutta法。具體求解過程如下。

由于式(20)為二階非線性微分方程，因此需要先將方程轉(zhuǎn)換為如下形式：

通過上式可將式(20)整理為

令如下等式成立：

則使用四階Runge-Kutta迭代求解上述方程，有

其中，δ為迭代時間間隔（δ=10-3）。這里采用python實現(xiàn)四階Runge-Kutta法，完成上述方程的求解。

2 數(shù)值模擬方法

2.1 物理場模型

為得到毛細界面的可視化演化過程，采用數(shù)值模擬方法研究了毛細管內(nèi)界面的振蕩運動過程。采用Comsol軟件中層流兩相流模型將流動物理場與相場耦合，假設(shè)毛細運動為不可壓縮且互不混溶的牛頓流體，其動態(tài)模型使用Navier-Stokes方程。

質(zhì)量守恒方程為

動量守恒方程為

其中，u，p，t分別表示速度場、流場中的壓力和時間；界面處的密度ρ和黏度μ分別根據(jù)兩相流體的體積分數(shù)加權(quán)計算得到，有

其中，Vf,1為流體1的體積分數(shù)，ρ1和ρ2分別為流體1和流體2的密度，μ1和μ2分別為流體1和流體2的黏度。

相場法是一種擴散界面方法，用于描述兩相流界面的動態(tài)變化過程，可以通過求解相場方程直接捕捉運動界面。相場方程為

其中，?，u，γ，λ，ψ分別為連續(xù)相場變量、流體速度、移動參數(shù)、能量混合密度、相場變量；?pf為控制界面厚度的參數(shù)。

移動參數(shù)γ為一個標量，定義為

其中，χ為遷移率，這里取χ=1m·s·kg-1。

混合能量密度λ有助于根據(jù)模型界面的表面張力調(diào)整界面厚度，定義為

相場變量ψ定義為

其中，?=1為其中一相，?=-1為另一相。

2.2 幾何模型及計算方法

建立的二維毛細管幾何模型如圖3所示。其中下部儲液罐為內(nèi)徑120 mm、高度45 mm的圓柱形容器，內(nèi)部盛有高度為30 mm的液體。計算中采用三種不同的流體工質(zhì)，具體參數(shù)見表1。數(shù)值模擬中采用的毛細管幾何結(jié)構(gòu)長度z= 100 mm，半徑R=1 mm，其中一端浸入液體中（z0=10R），另一端與大氣相連。圖3中藍色邊界設(shè)置為潤濕壁，將其設(shè)置為Navier滑移邊界條件。

表1 流體參數(shù)屬性Table 1 Properties of the fluid parameter

圖3 傾斜毛細管幾何模型Fig.3 Geometric model of tilted capillary tube in a liquid bath

3 結(jié)果與分析

3.1 模型驗證

為驗證上述理論模型及數(shù)值模擬結(jié)果的正確性，首先將常重力下垂直毛細管內(nèi)的理論求解結(jié)果、數(shù)值模擬結(jié)果與文獻[5]中的實驗數(shù)據(jù)進行對比。工作流體使用HFE-7100，毛細管半徑R=1mm，毛細管浸入液體的高度z0=6R，傾斜角度α=0，即毛細管垂直于液面放置，結(jié)果如圖4所示。可以看出，重力環(huán)境會誘導(dǎo)毛細管內(nèi)界面產(chǎn)生振蕩，通過實驗、理論模型和數(shù)值模擬得到的界面運動趨勢基本相同。三種結(jié)果相比，實驗獲得的振蕩過程阻尼最大，振蕩衰弱過程最快，同時上下往復(fù)振蕩的非對稱性較為明顯。這是由于實際物理過程中，毛細管內(nèi)壁面的缺陷和非理想清潔狀態(tài)均會影響動態(tài)潤濕過程，并引起能量損耗。理論結(jié)果中由于考慮了毛細管端口液體吸入和流出的能量損失（見圖2），較好地顯示出界面振蕩周期和非均勻振蕩效應(yīng)。數(shù)值模擬結(jié)果與實驗出現(xiàn)一定偏差，這是由于模擬條件完全理想化，并且對入口處的局部能量損失考慮不足。整體而言，理論分析和數(shù)值模擬結(jié)果均能夠復(fù)現(xiàn)由重力引起的毛細界面瞬態(tài)振蕩過程，能夠反映振蕩幅值、峰谷位置、非對稱性等特征。下文主要基于理論模型和數(shù)值模擬研究分析傾斜毛細管及不同重力環(huán)境下的界面振蕩動力學(xué)行為。

圖4 垂直毛細管界面振蕩過程的實驗結(jié)果及理論分析與數(shù)值模擬結(jié)果對比Fig.4 Comparison of results of theoretical modeling,numerical simulation and experiment

3.2 不同重力水平對毛細管內(nèi)界面運動的影響

為研究不同重力水平對毛細管內(nèi)界面振蕩行為的影響，將半徑R=1mm的毛細管垂直液面放置（即α=0），浸入液面高度z0=10R，分別將重力水平設(shè)置為1g/6，3g/8，g和 2g。理論模型分析及數(shù)值模擬計算得到的瞬態(tài)界面過程如圖5所示。可以看出，毛細管內(nèi)界面的穩(wěn)定高度隨重力加速度的減小而增高。毛細管內(nèi)的界面上升主要受毛細力驅(qū)動，當(dāng)重力水平減小后，更多流體能夠吸入毛細管內(nèi)，界面達到新的高度后，流體重力與毛細力相平衡。同時，從曲線的瞬態(tài)發(fā)展過程可以看出，毛細管內(nèi)界面的運動速度隨重力作用的減小而增加，這是由于垂直方向上界面的受力不平衡增加了液柱上升的瞬時加速度。

圖5 不同重力下毛細管內(nèi)界面的運動情況（實線為理論結(jié)果，虛線為數(shù)值模擬結(jié)果）Fig.5 Meniscus movement under different gravity levels (Solid line is for theoretical result while the dashed line is for simulation)

由圖5還可以看出，毛細管內(nèi)界面的振蕩效應(yīng)隨重力水平的增大而增強。對于同一流體工質(zhì)，Oh=μ/(ρRσ)1/2不變，而Bo=ρgR2/σ（其中g(shù)為重力加速度）隨重力水平的增加而增大。當(dāng)Oh?Bo時，毛細現(xiàn)象受重力作用主導(dǎo)，毛細管內(nèi)界面的運動從初始的慣性狀態(tài)過渡到振蕩狀態(tài)[19]。因此，在相同條件下，重力作用越大，毛細管內(nèi)界面振蕩效應(yīng)越顯著。同時，毛細管內(nèi)的界面振蕩行為受工質(zhì)物性的影響，為此定義了界面振蕩率，即

其中，hmax為振蕩的最大高度，hmin為振蕩的最小高度，S為毛細界面最終的穩(wěn)定高度。 通過本文分析可知，在常重力環(huán)境中，HFE-7100的振蕩率為0.96，而去離子水與乙醇的振蕩率相近，均為0.49。三種流體的振蕩率與其對應(yīng)的Oh/Bo關(guān)系一致，這表明毛細管內(nèi)界面的振蕩行為與Oh/Bo有關(guān)，隨著Oh/Bo的增大，毛細振蕩效應(yīng)逐步弱化。

為進一步分析界面振蕩過程中的能量耗散機制，針對工作流體HFE-7100，分析其在不同重力環(huán)境下主要的能量耗散隨時間的演變過程，如圖6所示。從圖6可以看出，在毛細運動過程中由重力勢能引發(fā)的耗散（）占主導(dǎo)地位（占總能量耗散的80%以上），幾乎比其他能量耗散高1個數(shù)量級。所有的能量耗散均隨著重力水平的減小而增加。這是由于能量耗散與毛細管內(nèi)界面運動的速度直接相關(guān)。如前所述，重力水平越小導(dǎo)致界面運動速度越大，因此對應(yīng)的能量耗散越大。

圖6 HFE-7100在毛細管內(nèi)界面運動的能量耗散過程Fig.6 Energy dissipation during capillary oscillation for HFE-7100

3.3 傾斜角度對彎液面運動的影響

進一步研究傾斜角度對毛細管內(nèi)界面振蕩行為的影響。在常重力環(huán)境下考慮毛細管相對垂直方向傾斜角度α為30°，45°，60°三種情況，結(jié)果如圖7所示。從圖7可以看出，無論傾斜角度如何，毛細管內(nèi)界面最終穩(wěn)定的垂直高度均相等。然而毛細管內(nèi)界面的振蕩行為與傾斜角度有較大關(guān)聯(lián)，毛細振蕩效應(yīng)隨傾斜角度的增加而逐步弱化。以HFE-7100為例，當(dāng)毛細管傾斜角度為45°時，振蕩率為1.21；當(dāng)毛細管傾斜角度增加到60°時，振蕩率增加到1.32。這是由于毛細管傾斜角度越小，沿毛細管方向的重力分量對毛細管內(nèi)界面爬升的影響越大，振蕩效應(yīng)增強。由于物性差異，三種不同工質(zhì)表現(xiàn)出不同的振蕩特征。其中，HFE-7100密度較大、表面張力小，與重力效應(yīng)相關(guān)的界面振蕩較為顯著，但界面平衡高度小；乙醇有較大的表面張力，因此對抗界面振蕩的能力最強。

圖7 不同傾斜角度下彎液面運動情況Fig.7 Meniscus movement under different tilt angles

圖8為典型的界面振蕩演化過程。此工況采用HFE-7100流體，毛細管傾斜角度分別為α=45°和α=60°。由圖8可知，t=0.15s 時毛細管內(nèi)界面處于向上爬升階段，t=0.25s時毛細管內(nèi)界面處于下降階段，t=1s時毛細管內(nèi)界面處于穩(wěn)定階段。從圖8可見，無論處于哪個階段，毛細管內(nèi)界面均垂直于毛細管，但不完全對稱，這與重力效應(yīng)相關(guān)。同時，毛細管外也有顯著的毛細效應(yīng)，且右液面高度明顯高于左液面，這與文獻[20]研究的實心圓柱體傾斜浸入液面的現(xiàn)象一致。

圖8 HFE-7100在毛細管內(nèi)運動的瞬態(tài)過程Fig.8 Meniscus movement of HFE-7100

4 結(jié)論

針對變重力及毛細管傾斜角度對毛細管內(nèi)界面振蕩行為的影響過程開展研究，建立了毛細管內(nèi)界面瞬時上升高度的理論模型，并開展了數(shù)值模擬研究，通過與實驗結(jié)果對比驗證了理論模型和數(shù)值模擬的可靠性。研究結(jié)果表明：重力水平減小，毛細界面垂直方向上的不平衡力差變大，毛細管內(nèi)界面穩(wěn)定高度增高，運動速度增加，能量耗散增加；當(dāng)重力增大時，毛細管內(nèi)界面的振蕩行為增強。同時，在毛細運動過程中由重力勢能引發(fā)的耗散占主導(dǎo)地位，占總能量耗散80%以上。通過對比三種不同流體工質(zhì)，發(fā)現(xiàn)HFE-7100由于密度最大，其振蕩行為最顯著，振蕩率為0.96，而去離子水與乙醇的振蕩率接近，約為0.49。當(dāng)毛細管傾斜時，毛細管內(nèi)界面的振蕩行為由工質(zhì)物性和沿毛細管方向的重力分量共同決定，隨傾斜角度的增大而增強。