疏水表面振動液滴模態演化與流場結構的數值模擬*

葉欣 單彥廣

(上海理工大學能源與動力工程學院, 上海 200093)

為了研究疏水表面垂直振動液滴的運動規律, 本文建立了振動液滴的三維模型, 考慮了振動液滴的動態接觸角變化過程, 通過流體體積函數和連續表面張力(volume of fluid-continue surface force, VOF-CSF)方法實現了液滴受迫振動的數值模擬, 得到了液滴的四種模態(2, 4, 6和8)動態演化過程、內部流場結構以及動態接觸角的變化規律. 隨著振動加速度的改變, 液滴可表現出豐富的模態, 而具體模態依賴于振動加速度的頻率變化. 以此為基礎, 本文對液滴的內部流場結構做了進一步的分析. 在模態2和模態4時, 液滴內部流動從底部向上產生“Y”型流動, 而在模態6和模態8時呈現對稱的渦流動. 且共振模態階數越高, 液滴內部速度平均值越大. 液滴振動時的動態接觸角明顯偏離靜態接觸角, 表明液滴振動模型有必要考慮動態接觸角. 模擬結果與文獻中的實驗結果做了對比, 結果符合良好.

1 引 言

固著表面振動液滴在強化傳熱、調控蒸發自組裝過程以及微流控芯片等領域中有著重要的應用[1-7]. 早在19世紀就有學者開始研究液滴的振動, Kelvin[8]和Rayleigh[9]對共振頻率下引起無黏自由液滴的動態行為進行了研究. 1932年Lamb[10]提出了自由液滴在不同的共振模態下的振動頻率計算公式, 稱為Rayleigh共振頻率公式. Strani和Sabetta[11]在不計重力的情況下, 對附著于固體表面上的非黏性液滴的軸對稱振動進行了研究, 與Lamb的研究不同在于, Strani和Sabetta從理論上分析了釘扎在固體表面上的液滴振動, 發現釘扎在固體表面上的液滴的第一模態形狀類似于自由液滴的第二模態形狀; 研究還發現第n振動模態中的無量綱固有頻率可以表示為接觸半徑與液滴半徑之比的函數. 之后, 大量論文討論自由液滴振動的性質, 主要探討液滴的振動方向、共振頻率、運動方向等方面的特征[12-17].

為深入研究液滴的振動特性, 大多數文章通過實驗研究液滴的振動特性. Brunet等[14], Noblin等[15]及Dong等[16]對在傾斜表面、垂直及水平方向上振動的液滴運動特征進行了研究. 周建臣等[17]實驗研究了液滴在不同頻率下的垂直振動特性, 研究發現液滴對外界驅動的不同響應與接觸線的振蕩行為、變形程度密切相關. 近些年諸多學者通過實驗研究自由液滴表面出現重力波的振動頻率范圍[18-21]. Noblin等[18]研究了接觸角和滯后作用對液滴振動模態的影響, 實驗通過改變垂直振動的頻率和振幅觀察到接觸線固著-移動兩種類型的振動模態. Shin和Lim[19], Kim和Lim[20]以及Park等[21]針對疏水玻璃片上垂直振動液滴內部的流動規律進行了研究, 實驗研究得到液滴垂直振動的2, 4,6, 8四種振動模態, 在模態2和模態4下液滴內部呈現“Y”型流動, 而在模態6和模態8頻率下呈現對稱的渦流動, 此外通過粒子圖像測速法測量(particle image velocimetry, PIV)發現模態4的流速大于模態2, 模態6和8的流速幾乎相同.Ramos[22]針對超疏水表面微液滴在不同振幅下的振動特性, 得出超疏水表面微液滴可看作自由液滴的結論. 只有少數文章通過數值模擬研究固體表面上液滴的振動問題[23-25]. James等[25]實現了固體表面上單個液滴在低頻和小振幅垂直振動下的數值模擬, 研究發現振動導致自由液滴表面形成重力波, 當振幅高于臨界值時會引發液滴破裂.

基于以上實驗研究發現, 液滴內部流動具有清晰的三維特征, 并且接觸角變化是關鍵參數之一.目前關于液滴振動的研究主要是通過實驗測定振動頻率、振幅和周期等因素對液滴內部流動產生的影響[19-21]. 而在實驗研究中, 由于透鏡效應, 導致無法準確測量到液滴自由界面附近邊緣層的流動特征. 因此, 本文在考慮動態接觸角的情況下, 建立了三維液滴振動模型, 實現了液滴在疏水表面上垂直振動的三維數值模擬, 獲得了液滴在不同頻率下垂直振動的四種模態, 預測了不同共振模態下液滴的內部流場結構, 有助于彌補實驗測量時由于存在透鏡效應而導致液滴自由界面附近邊緣層流場不準確的局限性.

2 數學模型

2.1 物理問題及假定

當基底受到強迫振動時, 液滴隨基底垂直上下振動, 當驅動加速度超過某一臨界值時, 液滴的內表面和外表面由于振動的影響產生壓力差, 在壓力差的作用下形成表面波. 隨著振動的進行, 表面波疊加形成葉瓣, 導致液滴發生形變. 液滴在上下振動時, 其內部的流動受到幾何形狀的限制, 在重力和黏性力的共同作用下, 形成渦流動或回流. 故在不同頻率下的垂直振動, 液滴內部形成不同的流場結構.

根據已有的實驗研究[19-21]可知液滴在特定的頻率下會呈現共振模態. 當液滴處于共振模態時,表面波疊加形成葉瓣, 葉瓣與葉瓣之間會出現節點, 液滴的振動模態階數等于節點數. 在共振模態時液滴表面形成的葉瓣是以中心線為軸對稱分布,葉瓣大小和數量相同, 所以液滴內部的微流動是軸對稱的. 本文主要研究液滴在垂直振動下的2, 4,6和8四種共振模態及其內部微流動, 故對三維液滴的振動過程做出如下假設:

1)液滴內部的流動是不可壓縮的, 并且在垂直振動下是軸對稱的;

2)由于液滴的振動周期遠遠小于同質量同體積的液滴蒸發的時間, 所以液滴振動過程不考慮自然蒸發;

3)由于本次模擬的液滴半徑約為1.14 mm,小于毛細長度約為2.73 mm, 其中g是重力加速度,ρl是液滴密度,σ是表面張力),相較于表面張力, 重力對液滴形狀的影響可忽略不計, 因此將疏水表面上的微液滴視為球冠形.

基于以上假設, 本文采用VOF-CSF模型[26,27]對三維液滴的振動過程進行模擬.

2.2 VOF模型

本文模擬的液滴振動涉及氣液兩相和相界面的形變, 需要得到氣液兩相的交界面, 因此可以采用VOF模型[26]. 設定空氣為主相, 液態水為次相.氣液兩相共用一套動量方程, 通過引進相體積分數這一變量, 實現對每一個計算單元相界面的追蹤,氣液界面采用幾何重建法處理. 模型控制方程如下.

連續性方程:

其中αs為次相體積分數;vs為次相速度. 主相體積分數不能通過該方程求解, 主相與次相體積分數之和為1,

動量方程:

式中,p為靜態壓力;v是流體速度,F是表面張力所導致的源項;ρ和μ分 別代 表平 均 密度 和動 力黏度:

(4)式和(5)式中下標 p 和 s 分別表示主相和次相.

2.3 連續表面張力(CSF)模型

動量方程(3)式中的表面張力F采用Brackbill等[27]提出的CSF模型進行處理. CSF模型主要是將VOF計算中附加的表面張力處理為一體積力,把該體積力作為源項加入到動量方程中. 本文模擬中只涉及水和空氣相, 所以(3)式中的體積力F在CSF中的表達方程式如下:

式中,σ為氣液兩相表面張力;κs為次相與主相之間的界面曲率, 通過垂直于相界面的表面局部梯度計算得到:

式中,n為相界面單位向量,n通過次相體積分數梯度計算得到, 如下式所示:

在VOF模型中采用連續表面張力模型時需要給定壁面接觸角, 用于調整壁面附近單元表面的法向,壁面附近的實際單元的表面法向量表示為:

式中,nw和tw分別是壁面的單位法向量和切向量;θ為接觸角.

2.4 接觸角模型

(9)式取決于接觸角θ, Chernova等[28]指出,有必要考慮接觸角的動態變化, 才能對液滴的內部流動進行正確的數值模擬. 根據已有的實驗研究表明, 動態接觸角與接觸線移動速度、液體物性、靜接觸角等因素相關. 本文液滴在疏水表面上振動,其接觸角遲滯現象相對較小, 因此本文忽略接觸角遲滯現象, 為了能夠同時處理前進接觸角和后退接觸角, 采用Kistler經驗公式[29]來計算動態接觸角:

式中,θ0為初始接觸角;θD為運動過程中隨速度變化的接觸角,Ca為毛細數;u為接觸線變化速度.

2.5 網格劃分與邊界條件

如圖1所示, 計算區域為正方體空間, 尺寸為3 mm × 3 mm × 3 mm, 整個計算空間可滿足液滴在振動條件下的振動幅度范圍和計算要求. 體積為5 μL 的液滴被布置在XY平面的中心位置上,XY平面是接觸角為 1 15°的疏水表面, 周圍是空氣,其他邊界均設置為壁面邊界. 網格劃分采用非結構化六面體網格, 為了選取合適的最小的計算單元,嘗試不同大小的網格尺寸, 分別用0.02, 0.03, 0.04和0.05 mm的網格尺寸對液滴進行了模擬, 模擬得到的液滴形態和內部流動情況偏差較小, 因此選用網格尺寸為0.04 mm, 既能達到模擬精度, 又可節約計算時間.

圖1 三維模擬計算區域(3 mm × 3 mm × 3 mm)Fig. 1. 3D computational domain (3 mm × 3 mm × 3 mm).

底面采用無滑移邊界條件, 即液滴和底面接觸處的流體速度等于底面的速度. 計算工質為空氣和水, 空氣為主相, 液態水為次相, 通過改變振動加速度函數的共振頻率, 研究不同共振頻率對液滴形態變化和內部流動規律的影響. 在基板上施加一個振動加速度, 使液滴做垂直振動:

其中a為振動加速度;f表示液滴的共振頻率;A為振幅, 模態2, 4, 6和8的振幅A分別為0.022, 0.046,0.013和0.0035 mm.

Lamb[10]總結出自由微液滴共振模態下的不同共振頻率計算公式, 稱為Rayleigh共振頻率, 表示為

其中f表示液滴的共振頻率;V,ρ和σ分別指液滴的體積、密度和表面張力;n是不小于2的整數, 表示模態數. 將液滴參數(ρ= 998 kg/m3,V= 5 μ L ,σ= 0.0728 N/m)代入(14)式, 計算得到模態2, 4,6和8的共振頻率. 見表1可知, 由Rayleigh方程計算所得的共振頻率與實驗所得的共振頻率較為接近.

表1 共振頻率的理論值和實驗值對比Table 1. Comparisons of theoretical and experimental results for resonance frequency of a 5 μ L droplet.

3 結果與討論

3.1 振動模態

圖2是體積為5 μ L 的液滴放置在水平的115°疏水表面上, 模擬的垂直振動下液滴的四種模態形貌與實驗結果[20]的對比圖. 圖2(a)是液滴在疏水表面上的共振模態實驗結果, 圖2(b)是液滴共振模態的模擬結果, 模擬所用共振頻率為表1中Rayleigh方程計算值. 由圖可知, 模擬所得的液滴模態與實驗結果符合較好, 這驗證了數學模型的可靠性. 根據前文所提, 因振動形成的波會在液滴表面疊加形成葉瓣, 且共振模態階數即為節點數.圖2中為了正確識別葉瓣, 葉瓣節點處以箭頭突出顯示. 如圖2所示液滴在四種共振模態中呈現出豐富的形態變化, 表面形成的葉瓣數目和節點數也有所不同. 圖2中表明葉瓣數目隨著共振模態階數的增加而增加, 葉瓣尺寸卻隨之減小. 因此, 當液滴處于模態8時, 由于葉瓣數目較多而葉瓣又較小導致識別困難. 結果表明, 模擬得到的液滴模態與實驗中的液滴模態吻合良好.

圖2 疏水表面上液滴共振模態(2, 4, 6和8)的模擬結果與實驗結果對比 (a)實驗結果[20]; (b) 模擬結果Fig. 2. Comparisons of the resonance modes (2, 4, 6, and 8) of a vibrating water drop on a hydrophobic surface: (a) Experimental results[20]; (b) simulation results.

3.2 液滴振動模態的動態演化過程

圖3是模擬的液滴分別在2, 4, 6和8模態一個周期內的形狀變化與實驗結果[21]的對比圖. 由圖3可知, 模擬所得到的液滴形態的變化與實驗結果符合較好. 隨著共振模態階數的增加, 液滴的振動周期越來越短, 約為12.4, 4.0, 2.0和1.2 ms. 這一現象與共振頻率以及振動幅度大小有關, 模態階數越高, 共振頻率越大, 振動周期越短. 前一個周期的最后相位與后一周期的第一相位重合. 如圖3所示, 液滴在模態6和8中形貌變化更豐富, 在模態2和4中則變化的相對簡單. 因為在高階模態中, 共振頻率較大, 液滴表面形成的波不斷疊加形成的葉瓣數更多, 導致模態6和8中的液滴呈現更豐富的形貌變化.

圖3 液滴在2, 4, 6和8共振模態下一個振動周期內形貌變化的實驗[21]與模擬結果對比Fig. 3. Comparisons between the experimental[21] and simulation results of the droplet shape evolution during one cycle of the vibration at resonance modes 2, 4, 6 and 8.

3.3 液滴振動過程的內部流場結構

目前關于液滴垂直振動的研究主要是通過實驗測定液滴振動周期、形狀、內部流動方式, 由于實驗測量時存在透鏡效應, 導致無法準確測量液滴自由界面附近邊緣層的流動特征. 為了進一步了解液滴垂直振動時的運動規律, 模擬三維液滴在共振模態2, 4, 6和8下的運動情況, 獲得液滴內部的流場結構.

圖4 是液滴在模態2, 4, 6和8下內部流場的三維視圖, 可以看出在每種模態下, 液滴內部的微流動都呈現軸對稱的渦流動. 液滴在模態2和4下, 液滴中心流動向上呈現“Y”型, 周圍伴隨著渦流動; 而在模態6和8時, 液滴中心流動向上到達液滴頂部, 周圍同樣伴隨著渦流動. 為了進一步了解液滴內部的流場結構, 在圖5中對液滴內部的二維剖面進行分析并與實驗結果對比. 圖5是2,4, 6和8模態下液滴內部的速度場模擬結果與實驗結果[20]對比, 左側是液滴內部流動的可視化實驗結果, 右側是模擬的液滴內部速度流線圖. 隨著模態階數的增加, 多樣化的流動模式得到了體現.

圖4 在2, 4, 6和8模態下液滴內部的三維流場圖Fig. 4. The three-dimensional flow field inside the droplet at modes 2, 4, 6, and 8.

圖5 液滴在2, 4, 6和8共振模態下內部流動結構的實驗[20]與模擬結果對比 (a) 模態2; (b) 模態4; (c) 模態6; (d) 模態8Fig. 5. Comparisons between the experimental[20]and simulation results of the internal flow structure of the droplet during the vibration at modes 2, 4, 6 and 8: (a) Mode 2; (b) mode 4; (c) mode 6; (d) mode 2.

圖5 (a)是液滴在模態2振動時內部的速度場模擬結果. 液滴在振動開始時, 底端的流體受力向上流動, 頂端流體則因為重力影響向下流動, 并在液滴上部形成對稱的回流. 隨著振動時間的增加,液滴的主體流動由中間向上, 經過液滴上側葉瓣疊加的節點處, 再沿著液滴表面回到三相接觸線上,同時在重力和黏性力的共同作用下, 在液滴中心軸兩側形成對稱的渦流動. 圖5(b)顯示了模態4的內部流場, 與模態2流動方式類似, 內部流動皆是自底部中心向上的循環流動, 最后返回液滴的中心. 模態4中的流型更長, 且產生的循環流動位置和大小與模態2不一樣, 在液滴上部會形成一個渦環. 故液滴在模態2和模態4中內部流動都呈現“Y”型流動.

圖5(c),(d)分別是模態6和模態8的內部流動情況. 如圖所示, 與模態2和模態4的流型不同,以液滴中心線為軸形成兩個大尺寸的渦流動, 并且在模態6時, 液滴邊緣處出現渦流動, 而在模態8時, 液滴邊緣處則沒有出現渦流動. 出現這種情況的原因是液滴在模態8時表面形成的葉瓣較多且較小, 幾何空間不足以誘導形成回流. 因此, 在模態6和模態8中液滴的內部流動循環變大形成對稱的渦流動而不是“Y”型流動.

為了進一步探究液滴的內部流動規律, 對液滴中心底部的垂直速度進行了計算. 在實驗測量時,由于存在透鏡效應, 通常是選取某一特定區域測量速度值, 無法準確測量液滴內部完整區域的速度值[20,21]. 為了與實驗測量速度值對比, 選取了同樣區域計算了液滴的平均垂直速度, 表2為液滴在四種共振模態下周期內的垂直速度平均值的實驗值[20]與模擬計算值. 從表2中可以看出, 模擬值與實驗值呈現同樣的趨勢, 在模態4時液滴速度平均值約為模態2速度值的兩倍; 而模態6和模態8的速度差值不大. 模擬結果與實驗結果數值上存在一定的差距. 從速度大小和量級上來看, 模擬的液滴速度應該小于或接近基板速度. 通過估算基板速度發現模擬得到的液滴局部區域平均速度值和基板速度值接近, 因此理論上模擬結果是準確的. 實驗和模擬數據之間的差異可能主要是實驗條件和測量誤差導致.

表2 液滴在模態2, 4, 6和8下中心底部的平均垂直速度的實驗[20]與模擬對比Table 2. Comparisons between the experimental[20] and simulation results of averaged vertical velocity at the central bottom region of the droplet in modes 2, 4, 6, and 8.

圖6為液滴完整區域在2, 4, 6和8模態下的平均速度隨時間的變化曲線圖. 從圖6中可以看出, 四種模態下的液滴平均速度皆呈現正弦波動變化趨勢, 且平均速度的峰值按模態階數的高低排序, 在模態8時的液滴平均速度峰值達到最大. 在模態2時由于振動頻率較低, 液滴內部速度變化趨于平緩, 而其他高頻模態下液滴速度波動較大.

圖6 在2, 4, 6和8共振模態下液滴內部平均速度隨時間的變化Fig. 6. The variations of the velocity with time at modes 2,4, 6, and 8.

3.4 液滴振動過程中動態接觸角的變化

液滴在隨基板垂直振動過程中, 由于接觸角的變化量不為零, 因此有必要考慮接觸角變化對液滴振動形態演變的影響. 如圖7所示, 考慮動態接觸角模擬得到的模態2液滴振動幅值隨時間變化曲線(Num.1)和考慮靜態接觸角模擬得到的模態2液滴振動幅值曲線(Num.2)具有相似的趨勢, 但在峰值上有一定的差異, 曲線(Num.1)和實驗測得的幅值變化曲線(Exp.)更為接近. 不考慮動態接觸角就會導致模擬的液滴高度偏低, 從而影響準確模擬液滴形態演化. 所以, 模擬中考慮了動態接觸角的影響.

圖7 在模態2時液滴振動幅值隨時間的變化. Num.1: 動態接觸角; Num.2: 靜態接觸角; Exp.: 實驗[20]Fig. 7. The variations of droplet vibration amplitude with time at mode 2. Num.1: dynamic wetting angle; Num.2:static contact angle; Exp.: experiment[20].

圖8 給出了在4種模態下液滴的動態接觸角隨時間的變化規律. 從圖8中可以看到, 在四種不同模態下, 液滴的動態接觸角都是以先增大到波峰再減小到波谷這樣的趨勢變化. 當施加振動加速度后, 液滴的接觸線受振動作用, 液滴從接觸線附近開始加速運動, 動態接觸角開始增大, 此時動態接觸角大于靜態接觸角為前進接觸角. 隨著振動持續液滴表面形成波瓣導致液滴變形, 液滴邊緣的壓力逐漸增大, 并且受到黏性阻力和毛細力的阻礙作用, 液滴的運動開始減速, 前進接觸角開始減小,造成接觸線停止前進運動, 此時動態接觸角接近于靜態接觸角115°. 由于液滴邊緣積累了一定的壓力, 在這個壓力的作用下, 液滴開始做回復運動,并造成接觸線的加速后退, 動態接觸角變小, 此時動態接觸角小于靜態接觸角115°為后退接觸角; 隨著回復運動的進行, 液滴邊緣壓力變小, 接觸線的后退開始減速, 后退接觸角增大直至接近靜態接觸角. 隨著液滴振動的進行, 接觸線的前進與后退運動重復之前的變化規律, 動態接觸角也隨之變化.

圖8 在2, 4, 6和8共振模態下液滴動態接觸角隨時間的變化Fig. 8. The variations of the dynamic contact angle with time at modes 2, 4, 6, and 8.

從圖8中還可以看到, 在前進接觸角階段, 接觸線的運動速度比較大, 因此動態接觸角偏離靜態接觸角比較多. 在本文的計算條件下, 最大接觸角偏離靜態接觸角接近6°. 綜合圖7和圖8結果,考慮動態接觸角對于準確模擬液滴形態演化尤為重要.

3.5 液滴振動過程中潤濕面積的變化

液滴在振動過程中會經歷鋪展和收縮兩個過程, 圖9為振動周期內液滴潤濕面積在2, 4, 6和8共振模態下的變化曲線, 潤濕面積皆呈現正弦波動變化趨勢, 液滴初始潤濕面積約為3.4 mm2. 從圖9中可以看出, 振動開始時, 液滴的潤濕面積逐漸增大, 在大約時液滴的潤濕面積達到最大值(波峰), 但是接觸線并沒有停止運動, 而是開始收縮, 潤濕面積開始縮小; 在大約時, 液滴的潤濕面積達到最小值(波谷), 此時接觸線速度接近于零, 再之后潤濕面積又開始增大重復之前的鋪展和收縮運動.

圖9 在2, 4, 6和8共振模態下液滴潤濕面積隨時間的變化Fig. 9. The variations of the wetting area with time at modes 2, 4, 6, and 8.

4 結 論

本文依據液滴共振頻率與模態的關系以及運動液滴的動態接觸角模型, 實現了在115°疏水表面上施加強制性振動時, 放置在疏水性表面上的液滴的不同共振模態(2, 4, 6和8)形狀和液滴內部的流動方式的模擬. 結果表明, 模擬所得115°疏水表面上液滴振動的模態與實驗結果符合良好. 主要得出以下結論.

1) 液滴在特定的頻率下發生共振并且呈現特殊的模態形狀, 較高的共振模態階數有更多的節點葉瓣. 隨著模態數的增加, 液滴表面形成的葉瓣數目增加, 尺寸逐漸減小. 在四種共振模態下, 模態階數越高, 共振頻率越大, 周期越短.

2) 在所有模態中, 液滴的內部流動方式呈現軸對稱的. 在模態2和模態4時, 內部流動從液滴底部中心開始, 經過液滴上側葉瓣疊加的節點處,呈現“Y”型流動方式; 在模態6和模態8中, 液滴內部流動呈現對稱的渦流動方式. 且共振模態階數越高, 液滴內部速度平均值越大.

3) 動態接觸角對液滴振動幅值變化有一定的影響, 且動態接觸角明顯偏離靜態接觸角, 表明考慮動態接觸角對取得正確數值模擬結果的必要性.