楔形體上復合液滴潤濕鋪展行為的格子Boltzmann 方法研究*

張曉林 黃軍杰2)?

1) (重慶大學航空航天學院,重慶 400044)

2) (重慶大學,非均質材料力學重慶市重點實驗室,重慶 400044)

固壁上液滴的潤濕鋪展行為是自然界中普遍存在的現象,針對楔形體上的復合液滴,采用基于相場理論的格子Boltzmann 方法對其潤濕鋪展行為進行探究.通過理論分析和數值模擬,發現液滴潤濕面積隨接觸角、楔形體頂角的減小而增大,液滴也越容易分裂.處于理論分裂臨界狀態附近的液滴,在一定密度比、黏度比條件下將沿楔形體壁面分裂.基于模擬結果生成以密度比、黏度比為坐標的液滴分裂狀態相圖,比較發現相同條件下初始狀態為平衡態的復合液滴更不易發生分裂.另外模擬還表明非對稱界面張力及非對稱運動黏度比也是影響液滴分裂結果的重要因素.

1 引言

由于固液、液液和氣液分子之間的不同相互作用,固壁上的流動界面存在豐富獨特的自然現象,在固體表面的眾多性質中,浸潤性是一個非常重要的特性.探究固壁上流體的潤濕和鋪展,能促進人們對固液相互作用的理解,這對于相關基礎科學研究和前沿技術應用的發展都具有特殊意義.尤其是隨著仿生材料制造技術的進步,已經實現了獨特的液滴濕潤和鋪展行為,推動著自清潔表面[1]、涂裝[2,3]、微液滴定向操作[4]和油水分離[5]等相關領域的快速發展.

液滴在固體壁面上的潤濕鋪展行為研究,主要涉及兩類表面: 一類是光滑表面,如一般的平面和曲面;另一類是含奇異點的非光滑表面,如簡單的楔形/矩形槽面、臺面和復雜的微結構表面.對于一般光滑平面來說,對應的就是楊氏方程[6]所描述的理想情況,在忽略重力時液滴平衡呈球帽形狀.Sui 等[7]從理論和實驗方面分析流體與球形固體之間的潤濕結果,發現接觸角與楊氏接觸角一致,液滴或氣泡體積無窮大時,存在極限潤濕位置.Li 等[8]研究了圓錐狀基底的曲率比對其上液滴自發定向運動的影響,當曲率比超過臨界值時,液滴就會停滯.Han 等[9]研究了二維液滴在兩側壁接觸角不同的圓角V 形槽上的潤濕行為,發現液滴的形態除與槽體開角、接觸角相關外,還依賴于液滴體積和圓角半徑.Wang 等[10]采用格子Boltzmann 方法(LBM)對兩平行圓柱纖維上液滴的形態進行了模擬,確定了不同纖維間距和液滴體積下的完整形態圖,存在三種可能的平衡構型: 桶形液滴、液滴橋和液柱,并提出了液滴形態轉變的能量和力的解釋.對于非光滑的表面,Herminghaus 等[11]討論了液滴在楔形、矩形槽上的潤濕,發現同一結構表面可以存在多種液滴形態.Chang 等[12]從理論和實驗兩方面驗證了錐臺表面角點對液滴潤濕鋪展所起的抑制作用.Zhou 等[13]探究了液滴在狹窄矩形臺面上的潤濕鋪展.Ma 等[14]則考慮了液滴在圓形、三角形以及方形柱上的釘扎現象.對于其他微結構表面上液滴的潤濕鋪展研究[15?17],不再詳述.

上述均是固壁上單液滴的潤濕鋪展研究,而復合液滴的潤濕鋪展屬于三相流問題,相較于單液滴的二元流情形存在更多流動界面,同時流體和固體之間相互作用也更加復雜,因此相關數值研究更具挑戰性.Said 等[18]對平壁面上復合液滴的潤濕鋪展行為進行了實驗和數值探究.Weyer 等[19]研究了圓柱纖維上油水復合液滴的形態及液滴分離的相關影響因素.Zhang 等[20]分析了平壁面和毛細管中復合液滴的平衡構型.He 等[21]和Li 等[22]模擬了圓形界面上復合液滴的潤濕鋪展.Huang[23]研究了圓形和橢圓形界面上液滴的潤濕平衡形態,并討論了圓柱纖維上兩液滴的吞噬包裹過程.由于液滴的潤濕鋪展行為與固壁形態存在很強關聯,而目前還沒有楔形體上復合液滴潤濕鋪展行為的相關研究,因此我們將就其開展探究,以提高對固壁上液滴微操作(如復合微液滴的分裂)的認知.

2 數值方法及理論模型

LBM 在誕生的30 多年里已經取得了長足的發展,具有許多優點,諸如算法的簡單性、并行性,同時能夠處理復雜的邊界條件[24].由于LBM 的介觀特性,在處理多相流問題時展現出了相當的優勢.本文采用混合LB-有限差分方法[23]進行數值模擬研究,流場方程采用LBM 求解,基于Cahn-Hilliard 相場理論的界面演化方程,則采用二階有限差分法進行空間離散和二階Runge-Kutta 法進行時間推進.

2.1 數值方法

2.1.1 三元流體相場理論

對于三元流系統,自由能泛函定義為

其中系數ai為常量.σij為流體i和流體j之間的界面張力,γ1,γ2,γ3為與界面張力相關的參數:

系數ai,κi與參數γi、界面厚度W之間的關系為:ai=3γi/(4W) ,κi=3γiW/8 .化學勢可表示為

其中,

可進一步得到

其中δi=γT/γi.Cahn 和Hilliard 將擴散通量近似為與化學勢梯度成正比的量[25],得到Cahn-Hilliard方程(CHE).對于三元流體,有兩個獨立方程:

三相流體交界處的界面角φi(見補充材料A圖S1)滿足

在固壁上對于CHE,化學勢μi應滿足無通量邊界條件

nS為固壁上指向流體的單位法向量;而潤濕邊界條件用于求固壁內與流體點相鄰的虛擬點處的相序參數?,以計算??(x,t) 和?2?(x,t) .采用文獻[26]發展的幾何潤濕邊界條件,假設擴散界面附近?的等值線相互平行,同時接觸線處界面切線與接觸角一致,隨后建立與該切線平行的過虛擬點的特征線,并利用相序參數在界面法向滿足雙曲正切函數分布的特征來計算此點的?值.該方法在曲面附近無需采用復雜的插值計算,也不用判斷界面和網格的相對構型,實施起來較為方便.另外上述潤濕邊界條件使用的是文獻[20]中提出的加權接觸角:

其中θij(ij) 表示流體i-j界面與固壁之間的接觸角,處于流體i一側,具體見補充材料圖S1.接觸角θij和θji互為補角,即θij=π?θji.而且三個接觸角不相互獨立,由靜力學平衡滿足如下關系:

對于??的計算,基于LBM 中D2Q9 模型使用的各向同性離散格式,在楔形體邊界處采用二階精度:

在離開楔形體邊界的地方采用四階精度[27]:

而?2?的計算在全域均采用二階精度:

其中δx為網格尺寸,δt為時間步長,格子速度為c=δx/δt,el為格子速度矢量,wl為權系數.

2.1.2 LBM 求解Navier-Stokes 方程

使用單松弛碰撞模型時的格子Boltzmann 方程為[28]

通過相序參數可以計算流體密度和運動黏度:

其中ρi和νi分別為流體i的密度和運動黏度(相應的動力黏度ηi=ρiνi).宏觀物理量的計算為

通過Chapman-Enskog 分析,(13)式在宏觀尺度近似下列方程:

對于任意形狀的固壁邊界,采用插值回彈方法[29]求解未知的分布函數fl.為提高穩定性,模擬采用多松弛碰撞模型[30].

2.2 物理模型

本文主要研究楔形體上二維復合液滴(實際上對應第三維度無限長的復合液柱)的潤濕鋪展問題.在液滴尺寸較小、表面張力相對較大的情況下可以忽略重力的影響.物理模型如圖1所示,計算域為Lx×Ly的矩形,楔形體頂角為χ,兩側壁長b.初始狀態的復合液滴為Janus 狀液滴,由關于過楔形體頂點的豎直線對稱的兩個半圓形液滴組成,半徑為r0,其中左側液滴為流體1 (Fluid1),右側液滴為流體2 (Fluid2),周圍環境為流體3 (Fluid3).計算域四周均采用周期邊界條件.對稱參數條件下σ12:σ13:σ23=1:1:1,接觸角θ13=θ23(根據(9)式,θ12=90?),流體1 與流體2 的密度比=1,運動黏度比=1,流體1 與流體3 的密度比為,運動黏度比為,動力黏度比為.文中取ρr=ρ1,νr=ν1,Lr=r0,Ur=σ12/η1,則有Tr=r0η1/σ12,Re=,模擬中Re=100.部分非對稱情況的參數見4.4 節.

圖1 物理模型圖示Fig.1.Physical model illustration.

3 網格密度和計算域大小檢驗

模擬中將Lr用NL個網格離散,Tr用Nt個時間步離散,可得網格尺寸δx=Lr/NL、時間步長δt=Tr/Nt.相場模擬中有無量綱參數Cn=W/Lr=(W/δx)/NL,在給定NL的條件下,應保證界面厚度W足夠小,但W過小時又無法精確捕捉界面,同時考慮到計算資源的消耗,取W/δx=4 .為減小棋盤效應的影響[31],計算域采用奇數網格離散.為計算穩定,取P e=2×104.

3.1 網格密度檢驗

網格數量的選擇需要考慮計算資源消耗和計算精度兩方面因素,我們在接觸角θ13=105?,楔形體頂角χ=90?的情況下,考察網格密度對液滴潤濕鋪展模擬結果的影響.液滴潤濕平衡形態的理論解推導見補充材料A.對于靜態問題,計算域大小不影響液滴的最終平衡形態,選用Lx×Ly=5×5的計算域,網格密度為NL=30.2,40.2,60.2,80.2,100.2,相應的Nt=302,402,602,802,2004.不同網格密度下液滴平衡形態的模擬結果見圖2(對稱性圖中只給出一側液滴的情形,后同),可見除三相點和接觸線位置附近的差別稍大外,界面其余位置的結果基本一致.我們還監測了接觸線位置隨時間的變化過程,由于初始時刻接觸線處于楔形體頂點處,因此以接觸線到楔形體頂點的距離dCL來表征接觸線的位置,相應的計算結果如圖3(a)所示.在液滴未分裂時,dCL還反映了其在楔形體表面的潤濕面積.相較NL=30.2,40.2,在NL=60.2,80.2,100.2 時,接觸線位置隨時間的變化比較接近,而在趨近平衡時,5 種網格密度下液滴的接觸線位置基本一致,綜合考慮,在靜力學計算中采用NL=40.2.對于動力學問題,針對Lx×Ly=5×5的計算域考察了=10,=10的液滴的動態潤濕鋪展過程,此時還采用了更大的網格密度NL=120.2,相應的各種網格分別對應Nt=302,402,1204,1604,2004,2404.液滴接觸線位置的演化如圖3(b)所示,可以看出在NL=100.2,120.2時的計算結果很接近,總體趨于收斂.但采用較大的網格密度對計算資源消耗巨大,同時接觸線位置演化的誤差主要存在于部分時刻(dCL為極值附近),且NL=60.2,120.2 時的相對誤差約為2.1%,故后面動力學問題中的計算采用NL=60.2.

圖2 不同網格密度下液滴的平衡形態Fig.2.Equilibrium morphology of the droplet at different mesh densities.

圖3 不同網格密度下液滴接觸線位置的演化 (a) ==1;(b) ==10Fig.3.Evolution of the contact line position of droplet under the different mesh densities: (a) ==1 ;(b) ==10 .

3.2 計算域大小檢驗

對于液滴分裂動力學問題,計算域大小將會影響液滴分裂的模擬結果,因此我們關注=10和=10的液滴在χ=90?和θ13=80?的楔形體上的潤濕鋪展,考察了5 種計算域大小(Lx×Ly=5×5,7.5×5.0,10×5,12.5×5.0,10×10,對應的網格數量分別為Nx×Ny=301×301,451×301,601×301,751×301,601×601)對模擬結果的影響.圖4 為上述5 種計算域中模擬得到的液滴最終平衡形態,比較發現計算域為Lx×Ly=7.5×5.0,10×5,12.5×5.0,10×10 時的液滴形態一致,且較Lx×Ly=5×5 的結果更準確.圖5 為上述不同計算域尺寸下液滴潤濕鋪展過程中接觸線位置隨時間演化的曲線,可見對于液滴動態演化過程也有類似的結果.因此取計算域x軸方向尺寸Lx≥7.5,同時y軸方向尺寸Ly≥5 對動力學模擬是合適的,本文采用Lx×Ly=10×5 的計算域.

圖4 不同計算域 Lx×Ly 下液滴的平衡形態Fig.4.Equilibrium morphology of droplets in different computational domains of Lx×Ly .

圖5 不同計算域Lx×Ly 下液滴接觸線位置的演化Fig.5.Evolution of the contact line position of the droplet in different computational domains of Lx×Ly .

4 結果與分析

4.1 楔形體頂角、壁面潤濕性對液滴平衡形態的影響

本節關注液滴潤濕的平衡態問題,主要討論楔形體頂角及壁面潤濕性對未分裂復合液滴平衡形態的影響.因流體物性對本節考慮的液滴最終平衡形態無影響,為簡單起見,模擬中選用的流體密度比、運動黏度比均為1∶1∶1 (即=1 和=1),計算域大小為Lx×Ly=5×5,其他相關參數為NL=40.2,Nt=402,C n=0.1 .首先探究楔形體頂角一定時,不同接觸角下液滴平衡形態的差異.圖6 為楔形體頂角χ=90?時隨接觸角變化的液滴平衡形態圖,能看到此時液滴的潤濕特性與平壁面類似,即接觸角越大,壁面越疏水,液滴潤濕面積越小,三相點到楔形體頂點處的距離就越大.在此給出各種接觸角下液滴三相點到楔形體頂點處的距離hnum,并與對應的理論值htheo(見補充材料A)進行比較,發現接觸角θ13接近臨界值2π/3?χ/2和 2π/3+χ/2 時,計算誤差稍大,而其余情況則比較吻合(圖7).

圖6 不同接觸角下的液滴平衡形態Fig.6.Equilibrium morphology of droplets at different contact angles.

圖7 h 隨接觸角的變化曲線Fig.7.Variation of h with the contact angle.

需要注意的是,接觸角為理論臨界分裂值時(2π/3?χ/2=75?),得到hnum>0,此時液滴并沒有分裂.我們還討論了楔形體頂角對液滴平衡形態的影響,圖8 為θ13=90?時三種不同頂角 (χ=75?,90?,105?) 楔形體上液滴的平衡形態,對應的液滴潤濕面積s分別為1.82,1.72,1.66,可知在壁面潤濕特性一定時,楔形體頂角越小,被潤濕的固壁面積就越大.根據理論分析我們在各種楔形體頂角、接觸角下給出了液滴的臨界分裂界線(圖9),當液滴處于臨界分裂界線左下方時對應分裂狀態,處于臨界分裂界線右上方時對應不分裂狀態.(90?,80?)表示楔形體頂角χ=90?,接觸角θ13=80?,此時理論預測液滴未分裂.分析表明接觸角、楔形體頂角越小,液滴越容易處于分裂狀態.

圖8 不同頂角楔形體上液滴的平衡形態Fig.8.Equilibrium morphology of droplets on wedges with different vertex angles.

圖9 液滴分裂/不分裂理論臨界界線Fig.9.Theoretical critical boundary of droplet splitting/non-splitting state.

4.2 密度比、黏度比對液滴潤濕分裂行為的影響

從4.1 節的分析可知,接觸角和楔形體頂角是影響液滴潤濕行為的兩個重要因素,而對于動態潤濕鋪展過程,當接觸角處于理論臨界分裂值附近時,其動力學結果可能不同,故本節探究楔形體頂角χ=90?,接觸角θ13=80?下密度比、黏度比對復合液滴潤濕分裂行為的影響.計算域大小為Lx×Ly=10×5,取NL=60.2,Nt=1204,則C n=0.07,Nx×Ny=601×301.圖10 給出兩種典型的液滴不分裂和分裂的潤濕鋪展過程,分別對應=50,=1和=50 ,=5 .從圖10 可以看到在演化的初期,液滴界面發生了扭曲變形,當三相點接近楔形體頂點時,在=1 的情況下液滴發生了回縮,最終不能分裂,而=5 時的液滴則發生了分裂.從速度場分布圖10 中可知,在演化過程中液滴界面附近產生了渦,隨著演化趨于平衡,渦逐漸遠離液滴,并隨黏性耗散而衰弱,直至消失.考慮液滴分裂與不分裂兩種結果,在圖11(a)建立了不同密度比、運動黏度比下楔形體上復合液滴潤濕鋪展的分裂狀態相圖,發現在運動黏度比≥3時,考慮的大部分密度比下的液滴發生了分裂,只是在小密度比=1 時,液滴沒有分裂,分裂發生在運動黏度比稍大的情況下.可見當較小時,液滴慣性較小,不利于實現液滴分裂,同時液滴密度一定時,在較大黏度比的情況下,環境流體黏度較小,對液滴潤濕鋪展的阻礙作用就越小,則液滴更容易分裂.圖11(a)稍做調整便得到了密度比、動力黏度比條件下的液滴分裂狀態相圖,如圖11(b)所示.

圖10 初始狀態為非平衡態Janus 狀液滴的潤濕鋪展過程及速度場分布 (a) =50 ,=1 ;(b) =50 ,=5Fig.10.Wetting and spreading process and velocity field distribution of Janus-like droplet with non-equilibrium initial state:(a) =50,=1 ;(b) =50 ,=5 .

圖11 初始狀態為非平衡態Janus 狀液滴的潤濕分裂狀態相圖 (a) rv13 -rρ13;(b) rη13 -rρ13Fig.11.Split/Non-split phase diagram of Janus-like droplets with non-equilibrium initial state: (a) rv13 vs. rρ13;(b) rη13 vs. rρ13.

在上述探究參數范圍內,我們還考慮液滴潤濕演化過程中的能量變化.圖12 為不同密度比液滴在潤濕鋪展過程中所達到的系統動能最大值Ek,max隨運動黏度比的變化曲線(小圓圈代表分裂,“+”符號代表未分裂),其中插圖為演化過程中系統動能的變化曲線,其表達式為

圖12 初始狀態為非平衡態時系統動能最大值隨運動黏度比的變化曲線Fig.12.Variation of the maximum kinetic energy of the system with the kinematic viscosity ratio under the nonequilibrium initial state.

4.3 液滴初始形態對其潤濕分裂行為的影響

前面討論的是初始條件為非平衡態Janus 狀液滴的情形,本節將探究初始形態為平衡態的復合液滴在楔形體上的潤濕鋪展行為.平衡態的復合液滴由非平衡態Janus 狀液滴演化而來,其半徑為r=0.7884r0,具體推導見補充材料B.初始時刻平衡態復合液滴的下三相點位于楔形體頂點處,其余參數與4.2 節一致,我們仍關注該情況下密度比、黏度比對液滴潤濕分裂行為的影響.圖13 為此初始條件下液滴分裂、不分裂的兩種典型過程,分別對應.與圖10 比較,發現=5時初始形態為非平衡的液滴發生了分裂,而初始形態為平衡態的液滴卻沒有分裂,分裂發生在=10 .以密度比、黏度比為坐標軸生成的液滴分裂狀態相圖如圖14 所示,比較圖11(a)和圖14(a)可知,此時液滴分裂發生在運動黏度比更大的情況;比較圖11(b)和圖14(b),發現復合液滴的分裂界線相較于初始狀態為非平衡態的情況更加緩和,說明在相同密度比的條件下,初始狀態為平衡態的復合液滴的分裂發生在動力黏度比更大的情況.總之,對于初始形態為平衡態的液滴,其界面能較低,系統總勢能較小,在相同的條件下更不容易突破能量壁壘而發生分裂.要使液滴分裂,應該提高潤濕鋪展前期的慣性作用,可以通過減小環境流體(或增大液滴)的密度或黏度來實現.

圖13 初始狀態為平衡態復合液滴的潤濕鋪展過程及速度場分布 (a) rρ13=50 ,rν13=5 ;(b) rρ13=50 ,rν13=10Fig.13.Wetting and spreading process and velocity field distribution of compound droplet with equilibrium initial state:(a) r ρ13=50,r ν13=5 ;(b) r ρ13=50 ,r ν13=10 .

圖14 初始狀態為平衡態復合液滴的潤濕分裂狀態相圖 (a) rv13 -rρ13;(b) rη13 -rρ13Fig.14.Split/non-split phase diagram of compound droplets with equilibrium initial state: (a) rv13 vs. rρ13;(b) rη13 vs. rρ13.

4.4 非對稱條件下復合液滴的潤濕鋪展行為

前面討論的均是對稱參數下復合液滴的潤濕鋪展行為,本節研究非對稱參數下復合液滴的潤濕鋪展.在研究分裂問題時,根據4.3節的結果,將討論的參數選取在分裂狀態相圖分界線附近.考慮界面張力不對稱時的兩種情形,分別對應σ12:σ13:σ23=1.000:1.115:0.816(?1=135?,?2=105°,?3=120?)以及σ12:σ13:σ23=1.000:1.155:0.577 (?1=150?,?2=90?,?3=120?),平衡時兩側液滴具有形狀上的非對稱性,相關推導見補充材料B.在此我們采用稍小的接觸角θ13=θ23=75?,根據(9)式可得θ12分別為 85.564?和 81.406?.初始時刻平衡態復合液滴的下三相點位于楔形體頂點處,其余參數與4.2 節一致.在對稱和非對稱界面張力的條件下,我們觀察了rρ13=50 和rν13=7 液滴的潤濕鋪展過程,其左右側接觸線位置隨時間的演化如圖15 所示(圖例中LS 表示左側液滴接觸線,RS 表示右側液滴接觸線).在該條件下發生分裂的臨界值為dCL=2.35,當dCL>2.35 時說明液滴發生了分裂.圖15 還給出t=0,2,4,30 時刻液滴的形態,可以看到在σ12:σ13:σ23=1:1:1 和σ12:σ13:σ23=1.000:1.115:0.816兩種情況下,液滴均發生了分裂,但相較于對稱情形,非對稱界面張力下液滴的分裂程度較低;在σ12:σ13:σ23=1.000:1.155:0.577時,液滴沒有分裂.根據上述三種界面張力下的模擬結果,我們發現隨某一界面角的減小,一側液滴被另一側液滴包裹的程度更大,相應地就越不容易發生分裂.

圖15 左右側液滴接觸線位置的演化Fig.15.Evolution of the position of the left and right droplet contact lines.

圖16 左右側液滴接觸線位置的演化Fig.16.Evolution of the position of the left and right droplet contact lines.

5 結論

本文對楔形體上復合液滴的潤濕鋪展行為進行了數值研究.對于固壁上液滴潤濕的靜態問題,理論分析表明接觸角、楔形體頂角越小,液滴越容易分裂.在不發生分裂的情況下,楔形體頂角一定時,液滴的潤濕特性與平壁面類似,即接觸角越大,壁面越疏水,液滴潤濕面積越小;壁面潤濕特性一定時,楔形體頂角越大,被潤濕的固壁面積就越小.對液滴潤濕鋪展的動態問題,當液滴處于理論臨界分裂接觸角附近時,在一定的密度比、黏度比條件下會沿楔形體壁面發生分裂.基于模擬結果生成以密度比、運動黏度比為坐標的液滴分裂狀態相圖,通過能量分析表明密度比、運動黏度比越大時,潤濕鋪展過程中的慣性效應越大,液滴就更容易分裂.對比液滴的分裂狀態相圖,發現在相同的條件下初始形態為平衡態的復合液滴更不易發生分裂.對于非對稱情況,我們發現左右側液滴運動黏度差的增加有利于液滴分裂,同時界面張力不對稱導致的液滴被包裹程度越大,液滴越不易分裂.本研究只討論了對稱參數和部分非對稱參數下楔形體上二維復合液滴的潤濕鋪展行為,對于三維問題,將在以后的工作中進一步探討.