尹宗軍，張椿英，蘇 蓉，徐 輝，陳倩楠

(安徽信息工程學院 機械工程學院，安徽 蕪湖 241100)

0 引言

流體力學表明，液滴潤濕時流體和固體邊界之間存在滑移運動.滑移運動意味著固液界面處的速度不連續(xù)性，而滑移長度λ是一個假設的邊界內(nèi)距離，在該距離處，固液界面處的流體速度達到固體邊界速度.吳承偉等人[4]綜述了液體流動的邊界滑移及其相關問題的發(fā)展過程，探索了邊界滑移的產(chǎn)生機理以及各因素對邊界滑移的影響規(guī)律.近年來，相關文獻表明粗糙壁面、微通道、外部氣體環(huán)境、高壓高剪切以及驅(qū)動壓強差對邊界滑移長度都具有較大影響[5-15].趙士林[16]認為滑移現(xiàn)象對超疏水表面材料的研究具有重要意義.秦長劍[17]采用耗散粒子動力學的方法，研究了微通道的邊界條件及表面潤濕性對邊界滑移的影響.曹炳陽[18]研究了滑移現(xiàn)象對納米尺度液體流動的影響規(guī)律及其機制.許少鋒等人[19]研究了微通道疏水表面的滑移現(xiàn)象，模擬了平板間的Couette流動.他們的研究結果表明壁面與流體間排斥作用越強，即疏水性越強，壁面滑移越明顯，并且滑移長度與接觸角之間存在近似的二次函數(shù)關系.黃元丁[20]認為表面微觀形貌對固體表面的潤濕性和黏滯力有重要影響.他們著重探討了表面微觀形貌對邊界滑移的影響，研究了五種狀態(tài)的Cassie態(tài)表面上表面形貌、氣液接觸比、流道高度與驅(qū)動壓強差對邊界滑移長度的影響.

近年來，由于實驗研究成本過高，數(shù)值仿真技術被大量運用到解決復雜流體問題中.VOF、CLSVOF(Coupled Level Set and Volume of Fluid)、相場動力學等方法分別被運用到液滴撞擊領域中，模擬得到的結果與實驗取得高度一致性.值得一提的是，模擬手段對改變撞擊參數(shù)更為簡單，且能給出液滴內(nèi)部的動力學信息，例如壓力、速度等，這將有利于某些流動現(xiàn)象的機理解釋.然而，已有的研究多數(shù)是基于靜態(tài)網(wǎng)格，也就是網(wǎng)格化后，網(wǎng)格不再變化，這不利于捕捉像液滴這樣形變較大的流體.若要使用靜態(tài)網(wǎng)格精確捕捉，就需要細化整個計算域，這不利于減少計算時間.動態(tài)自適應網(wǎng)格可以精確捕捉所設定的物理量的變化，在物理量變化劇烈的地方自動細化網(wǎng)格，而在物理量變化不劇烈的地方自動粗化網(wǎng)格加快計算.因此，當下迫切需要解決的技術問題是：(1)如何使用模擬計算取代實驗仍然能獲得精確的結果；(2)如何在模擬計算中用最少的網(wǎng)格精確獲得界面形變，節(jié)約計算物理時間.采用VOF界面追蹤方法捕捉相界面，利用動態(tài)網(wǎng)格自適應手段減少不必要的計算，且在相界面處自動加密，就為單液滴撞擊平坦固體表面的模擬提供了可行的手段.因此，本文基于VOF界面追蹤和自適應動態(tài)網(wǎng)格技術探究了液滴撞擊壁面的振蕩行為，進一步揭示固液界面邊界滑移對液滴振蕩特性的影響.

1 控制方程與參數(shù)設置

1.1 控制方程

在VOF模型中，不同的流體組分共用一套相體積分數(shù)方程，通過引進相體積分數(shù)這一變量，實現(xiàn)對每一個計算單元相界面的追蹤.液滴撞擊表面模型只涉及液滴及空氣兩種流體，假設二者之間為不可壓縮、互不融合的層流運動，定義液相液滴為主相(ρL和μL用于表示液滴的密度和黏度，相分數(shù)設為c=1)，氣相空氣為第二相(ρV和μV分別表示空氣的物理特性，相分數(shù)設為c=0)，而c介于0和1之間時代表兩相界面區(qū)域；即可將兩相流動等效為一種混合物(密度ρ和黏度μ).該混合物的密度ρ和黏度μ可以用體積分數(shù)變量c分別表示為ρ=cρL+(1-c)ρV，μ=cμL+(1-c)μV.針對兩相流流體運動的特點，流動控制方程包括連續(xù)性方程，帶表面張力的動量方程以及VOF相分數(shù)方程[21]：

▽·u=0,

(1)

(2)

(3)

式中，u≡(u,v)為速度矢量，p為壓強，D定義為Dij≡(?iuj+?jui)/2，g為重力加速度，δ為迪拉克算子，κ表示界面的平均曲率，n表示從液相流出的界面的單位法線.

1.2 離散方程

▽·un=0,

(4)

(5)

(6)

上述離散方程可進一步采用古典時間分解投影算法簡化為[22]：

▽·un+1=0,

(7)

(8)

(9)

(10)

式中，ρi為密度，ui、μi、Di以及σi為對應時間節(jié)點i(i可取n、n+1/2和n+1)的值，u*表示輔助速度場.結合(7)式和(10)式，我們有以下Poisson方程：

調(diào)查結果顯示，有7%的學生參加了3項以上的課題研究，28%的學生參與了1～2項的課題研究。碩士研究生參與課題研究對其科研水平的提升有很大的幫助，經(jīng)常參與科研課題，能夠激發(fā)學生的學習積極性和科研創(chuàng)新能力。

(11)

GERRIS軟件是求解描述流體流動的偏微分方程的自由軟件程序，其源代碼是免費的.它具有以下主要特點：(1)求解含時不可壓縮變密度Euler、Navier-Stokes方程；(2)自適應網(wǎng)格細化,分辨率根據(jù)流體的特征動態(tài)調(diào)整；(3)時間空間二階離散精度；(4)可以靈活的附加源項；(5)支持MPI并行、動態(tài)負載平衡、并行脫機可視化；(6)精確的表面張力模型.基于GERRIS軟件中成熟的VOF界面追蹤技術以及自適應動態(tài)網(wǎng)格手段完全可以實現(xiàn)對單液滴撞擊平坦固體表面的仿真模擬.(11)式使用GERRIS中的基于四元/八叉樹的多級解算器求解該方程.同時，(8)式可以重新改寫為

(12)

這是一個Helmholtz型方程，可以使用GERRIS中的多級Poisson解算器求解.

1.3 邊界條件

(13)

滑移邊界條件和無滑移邊界條件，這兩種邊界條件的法向流速都為零(即v=0)，區(qū)別在于這兩種邊界切向速度u是否為零；即無滑移條件就是指在靠近壁面處流體的流動速度為零u=0，而有滑移條件就是它們之間有相對速度.滑移邊界條件可以表示為

u=u0+λ?u/?y，

(14)

1.4 計算域與初始條件

首先在GERRIS軟件中定義理想的計算域尺寸，液滴撞擊壁面過程可以簡化為一個Lx×Ly=2L×L二維計算域，L=8 mm，如圖1(a)所示.計算域的右側、左側以及上側邊界指定為壓力出口邊界，底部指定為滑移邊界條件.在計算域內(nèi)部包含一個圓形區(qū)域(液滴，直徑D0)，液滴距離初始壁滴距離h=D0，以便液滴周圍的氣體流動能夠在液滴撞擊表面之前以物理方式發(fā)展，即保證將插入動量中的計算表面張力項始終具有最高的精確度.

液滴在鋪展過程中液滴的形態(tài)是時變的，這就要求網(wǎng)格分辨率沿液滴/氣體界面變化，因此動態(tài)自適應笛卡爾網(wǎng)格細化技術所用網(wǎng)格細化的判據(jù)可以表示為

‖c‖Δx/max‖u‖>δ,

(15)

式中，‖c‖是VOF體積分數(shù)變量c的值，Δx為子網(wǎng)格的尺寸，max‖u‖為流體穿過局部單元的最高速度，δ為閾值參量(δ取0.01).若(15)式滿足，則GERRIS軟件自動將網(wǎng)格細化一級.類似地，若滿足下式

‖c‖Δx/max‖u‖<δ/4,

(16)

GERRIS軟件自動將網(wǎng)格粗化一級，以減少計算.加密等級是指GERRIS軟件中計算域是由一個個子正方形格子組成，若給定加密等級為10是指對構成計算的子正方形格子進行210=1024次切分，也就是說如果正方形格子邊長為1 mm，則最小的網(wǎng)格長度為0.000 98.顯然，加密等級越高，網(wǎng)格越密，計算精度越高，但物理計算時間也就越長.液滴界面局部區(qū)域內(nèi)加密等級選擇10；其余區(qū)域加密等級選擇6，故相應的計算單元尺寸為L/210和L/26(即7.8125×10-3mm和1.25×10-1mm).初始的網(wǎng)格使用GERRIS軟件生成.

初始時，壓強p均設置為0，重力加速度g向下.如圖1(b)所示，應用的數(shù)值網(wǎng)格的放大區(qū)域顯示了移動的液滴以及網(wǎng)格的高分辨率，有助于更好地理解動態(tài)網(wǎng)格細化技術.圖1(c)給出了液滴撞擊過程中(i)碰撞、(ii)鋪展以及(iii)反彈三種運動特征，表明網(wǎng)格能自適應地跟隨液滴.

2 結果與分析

2.1 液滴撞擊光滑固體表面的運動特征

我們將水滴(D0=2.0 mm，V0=0.31 m/s)與光滑固體超疏水表面(靜態(tài)接觸角θE=157°)的VOF模擬與文獻[23]中的實驗進行了比較.在計算中，水和空氣的密度和黏度分別設定為ρL=9.98×102kg·m-3，ρV=1.2 kg·m-3；μL=1.003×10-3N·s·m-3，μV=1.8×10-5N·s·m-3，表面張力設置為γLV=7.5×10-2N·m-1.邊界條件為自由滑移邊界.

圖1 液滴的滑移邊界條件及運動特征

圖2 液滴撞擊光滑固體表面的運動演化

圖2(a)和(b)所示分別為VOF模擬與文獻[23]的液滴撞擊光滑固體表面的運動演化圖，圖2(c)顯示了鋪展因子β的時間變化曲線(定義為β=L/D0，L為當前的潤濕長度).當液滴接觸壁面時，液滴底端流體開始向四周鋪展.定義時間節(jié)點t=0.0 ms為水滴剛接觸到固體表面時刻，此時β=0.大約1.6 ms后，液滴底部出現(xiàn)褶皺型，由于疏水表面的黏附力影響，液滴的下部并沒有形成向四周的徑向流動.隨著液滴高度的降低，在t=3.8 ms時，液滴底部達到最大潤濕面積，最大鋪展因子βmax=1.85.在時間t=4.0 ms時，液滴呈現(xiàn)扁平狀，液滴中出現(xiàn)向上的突出.這種行為的出現(xiàn)實際上與能量轉(zhuǎn)換有關.當液滴撞擊到固體基底上時，由于固體表面的阻擋，液滴的法向動量發(fā)生變化，液滴周向鋪展以擴大表面積，直到動能完全轉(zhuǎn)化為表面能量；這種不平衡狀態(tài)是暫時的，在表面張力的作用下，液滴已經(jīng)周向鋪展的流體會向中心回縮，形成反向慣性力.在時間6.8 ms到10.8 ms內(nèi)，這種反向慣性力迫使?jié)竦蚊娣e減小，β值不斷減小.在t=11.2 ms時，這個反向慣性力迫使水滴從基板上反彈.圖2(c)進一步表明液滴撞擊光滑固體表面過程中模擬和實驗的液滴鋪展因子β變化具有一致性.

2.2 固液界面邊界滑移對液滴振蕩的影響

為了使我們的結論更具普遍性，我們使用了兩組(V0=0.2 m/s,θE=70°和V0=0.4 m/s,θE=110°)來充分解釋滑移長度λ對液滴振蕩特性的影響.更具體地說，在VOF計算中我們使用了五種不同的邊界條件：λ=0.005 mm、0.2 mm、0.4 mm、0.8 mm以及自由滑移邊界.這里，選擇直徑D0=1.6 mm的水滴作為仿真模擬對象.

圖3 不同滑移邊界條件下鋪展因子β的時間演化

圖3顯示了在不同滑移邊界條件下，當V0=0.2 m/s,θE=70°(親水表面，近毛細管鋪展)和V0=0.4 m/s,θE=110°(疏水表面，近反彈)時，液滴在光滑表面上鋪展因子β的時間演化圖.圖3左側用于表示第一個振蕩周期(i)，這有利于描述液滴的最大鋪展因子βmax；圖3右側用于表示100 ms內(nèi)液滴的振蕩行為(ii)，這有利于描述液滴振蕩行為的基本物理量，如平均振幅比χ以及平均振蕩周期τ.在圖3中，無論是近毛細擴散還是近反彈，都觀察到了液滴振蕩變化的類似曲線.可以看出，近毛細鋪展比近反彈的振幅更大.圖3(a)-(i)和圖3(b)-(i)均表明，在第一個振蕩周期下，鋪展因子βmax的值略有變化.在第一個振蕩周期后，邊界滑移顯著改變了液滴振蕩的平均振幅比χ以及平均振蕩周期τ，尤其是在組V0=0.2 m/s,θE=70°中(見圖3(a)-(ii)).其原因是鋪展長度較大時，較小的滑移長度導致了更多的能量損失.圖3(a)-(ii)表明，當λ=0.05 mm變化到自由滑移邊界條件時，平均振幅比χ和平均振蕩周期τ分別為1.073/14.7 ms，1.066/15.7 ms，1.064/16.1 ms，1.056/16.7 ms和1.054/17.3 ms.λ=0.05 mm的χ值比自由滑移條件下的值小了2%，τ值減小了15%.這意味著隨著λ值的減小，振蕩階段消失得更快.對于圖3(b)-(ii)，平均振幅比χ和平均振蕩周期τ的值分別為1.074/13.7 ms，1.073/13.7 ms，1.066/14.1 ms，1.064/14.2 ms和1.063/14.2 ms.圖3表明，邊界滑移λ值越小，固液界面處的液體剪切速率越顯著，這將更有效地抑制液滴振蕩.

3 結論

當前的研究針對帶有移動接觸線的氣液兩相流動的數(shù)值研究，在固壁邊界處引入滑移邊界條件.通過GERRIS軟件內(nèi)部二階精度的離散格式，利用成熟的VOF界面追蹤技術以及自適應動態(tài)網(wǎng)格手段實現(xiàn)對液滴垂直撞擊光滑表面進行了二維VOF模擬.本研究的主要目的在于探討固液界面邊界滑移對液滴振蕩行為的影響.具體而言，是為了了解液滴撞擊固體表面的振蕩行為如何響應滑移邊界的變化，并提供一些可行的策略來縮短一些潛在應用的振蕩時間.研究結果表明在較低的滑移長度值下，如λ<0.05 mm，液滴在親水表面上鋪展，鋪展因子β衰減變得更加明顯.在滑移長度較高的情況下(包括自由滑移條件)，如λ>0.2 mm，滑動長度的增加對親水性或者疏水性表面上的鋪展特性均影響不大.