固液界面邊界滑移對(duì)液滴振蕩行為影響的VOF模擬

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

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

0 引言

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

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

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

1.1 控制方程

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

▽·u=0,

(1)

(2)

(3)

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

1.2 離散方程

▽·un=0,

(4)

(5)

(6)

上述離散方程可進(jìn)一步采用古典時(shí)間分解投影算法簡(jiǎn)化為[22]：

▽·un+1=0,

(7)

(8)

(9)

(10)

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

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

(11)

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

(12)

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

1.3 邊界條件

(13)

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

u=u0+λ?u/?y，

(14)

1.4 計(jì)算域與初始條件

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

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

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

(15)

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

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

(16)

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

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

2 結(jié)果與分析

2.1 液滴撞擊光滑固體表面的運(yùn)動(dòng)特征

我們將水滴(D0=2.0 mm，V0=0.31 m/s)與光滑固體超疏水表面(靜態(tài)接觸角θE=157°)的VOF模擬與文獻(xiàn)[23]中的實(shí)驗(yàn)進(jìn)行了比較.在計(jì)算中，水和空氣的密度和黏度分別設(shè)定為ρ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，表面張力設(shè)置為γLV=7.5×10-2N·m-1.邊界條件為自由滑移邊界.

圖1 液滴的滑移邊界條件及運(yùn)動(dòng)特征

圖2 液滴撞擊光滑固體表面的運(yùn)動(dòng)演化

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

2.2 固液界面邊界滑移對(duì)液滴振蕩的影響

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

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

圖3顯示了在不同滑移邊界條件下，當(dāng)V0=0.2 m/s,θE=70°(親水表面，近毛細(xì)管鋪展)和V0=0.4 m/s,θE=110°(疏水表面，近反彈)時(shí)，液滴在光滑表面上鋪展因子β的時(shí)間演化圖.圖3左側(cè)用于表示第一個(gè)振蕩周期(i)，這有利于描述液滴的最大鋪展因子βmax；圖3右側(cè)用于表示100 ms內(nèi)液滴的振蕩行為(ii)，這有利于描述液滴振蕩行為的基本物理量，如平均振幅比χ以及平均振蕩周期τ.在圖3中，無論是近毛細(xì)擴(kuò)散還是近反彈，都觀察到了液滴振蕩變化的類似曲線.可以看出，近毛細(xì)鋪展比近反彈的振幅更大.圖3(a)-(i)和圖3(b)-(i)均表明，在第一個(gè)振蕩周期下，鋪展因子βmax的值略有變化.在第一個(gè)振蕩周期后，邊界滑移顯著改變了液滴振蕩的平均振幅比χ以及平均振蕩周期τ，尤其是在組V0=0.2 m/s,θE=70°中(見圖3(a)-(ii)).其原因是鋪展長(zhǎng)度較大時(shí)，較小的滑移長(zhǎng)度導(dǎo)致了更多的能量損失.圖3(a)-(ii)表明，當(dāng)λ=0.05 mm變化到自由滑移邊界條件時(shí)，平均振幅比χ和平均振蕩周期τ分別為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%.這意味著隨著λ值的減小，振蕩階段消失得更快.對(duì)于圖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 結(jié)論

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