沈?qū)W峰 曹宇 王軍鋒 劉海龍

(江蘇大學(xué)能源與動力工程學(xué)院, 鎮(zhèn)江 212013)

(2019 年 11 月 2日收到; 2019 年 12 月 31日收到修改稿)

1 引言

液滴撞擊固體壁面現(xiàn)象廣泛存在于工農(nóng)業(yè)領(lǐng)域, 如內(nèi)燃機中油滴撞擊缸壁、噴霧冷卻、噴涂印染、噴墨打印、農(nóng)藥在作物表面沉積等. 液滴撞擊壁面后的行為包括沉積、邊緣高速濺射、冠狀濺射、回縮破碎、部分回彈及完全回彈等[1]. 液滴的物性參數(shù)、撞擊基面的理化性質(zhì)以及撞擊條件共同決定著液滴撞擊固體壁面后的動態(tài)行為過程[2,3].

在工業(yè)生產(chǎn)過程中, 流體常因添加了活性劑、高分子材料或顆粒而表現(xiàn)出不同的非牛頓特性, 包括黏彈性、剪切增稠、屈服應(yīng)力和剪切變稀等. 牛頓流體液滴撞擊固體壁面行為已被學(xué)者們廣泛探討[4,5]. 然而, 由于非牛頓流體流變學(xué)行為的復(fù)雜性, 目前有關(guān)非牛頓流體液滴撞擊固體壁面行為的研究尚少. Jung等[6]與Huh等[7]的研究表明在流體中添加極少量的聚合物材料, 在不明顯改變流體的密度, 表面張力和剪切黏度等物性的情況下, 可使流體表現(xiàn)出顯著的黏彈性, 抑制液滴撞擊固體壁面后的回彈行為. Luu 等[8]和 Sa?di等[9]通過實驗發(fā)現(xiàn)對于表現(xiàn)出屈服應(yīng)力特性的非牛頓流體, 液滴撞擊固體壁面后的鋪展和回縮行為都受到了抑制.Boyer等[10]指出, 剪切增稠液滴撞擊固體壁面后的鋪展行為與牛頓流體液滴表現(xiàn)出極大的不同, 在撞擊壁面后的初始鋪展階段液滴的變形會突然停止. 對剪切變稀流體液滴撞擊固體壁面行為的研究, 主要采用的是實驗方法, 關(guān)于剪切變稀特性對液滴撞擊壁面行為的影響機制尚存在爭議. German等[11]研究了冪律流體液滴撞擊固體壁面后的鋪展直徑與液膜厚度的變化, 結(jié)果表明剪切變稀特性影響著液滴撞擊固體壁面后的動力學(xué)行為, 但液滴稠度系數(shù)在鋪展過程中起主導(dǎo)作用. 然而An等[12,13]對剪切變稀流體液滴撞擊不同疏水性壁面的過程研究后指出, 在撞擊過程中動態(tài)變化的剪切黏度直接影響著液滴的鋪展. Andrade等[14]的研究同樣指出, 撞擊到固體壁面后, 液滴內(nèi)部動態(tài)變化的流場會改變剪切變稀流體液滴的剪切黏度, 直接影響著液滴的動力學(xué)行為. Laan等[15]研究了表現(xiàn)出剪切變稀特性的血液液滴撞擊壁面的行為, 他們指出相較于無窮剪切黏度, 剪切變稀特性對液滴行為的影響可忽略不計. 劉海龍等[16]將納米顆粒分散到環(huán)氧樹脂中配制了表現(xiàn)剪切變稀特性的納米流體,研究了納米顆粒的加入對液滴撞擊固體壁面后鋪展行為的影響, 他們指出納米顆粒的加入顯著抑制了液滴在壁面上的鋪展和回縮行為, 但通過實驗難以研究流變學(xué)參數(shù)對液滴撞擊的獨立影響.

雖然測量與分析技術(shù)的發(fā)展已經(jīng)可以幫助我們通過實驗從唯象的角度認知液滴撞擊固體壁面這一復(fù)雜的物理現(xiàn)象, 但是由于實驗條件的局限性, 對流場范圍內(nèi)的一些物理量如速度與壓力仍然難以精確測定, 如相界面移動等微觀物理現(xiàn)象也難以捕捉. 隨著計算機技術(shù)和數(shù)值計算方法的不斷進步, 數(shù)值模擬作為一種基本手段, 被眾多學(xué)者用來研究液滴撞擊固體壁面行為. 由于氣體和液體兩相間存在的高黏度與密度比, 處理撞擊液滴的相界面移動問題是采用數(shù)值模擬方法研究液滴撞擊固體壁面行為的難點. 在固定網(wǎng)格上求解相界面移動的方法主要分為兩類: 顯式求解(界面追蹤法)和隱式求解(界面捕捉法). 邊界積分法[17]和鋒面追蹤法[18]是兩種典型的界面追蹤方法, 界面追蹤方法提供的數(shù)值精度較高, 但它的適用范圍僅限于斯托克斯流動. 與界面追蹤方法不同, 界面捕捉法通過求解固定網(wǎng)格上獨立的相函數(shù)來獲得相界面的移動. 目前被學(xué)者應(yīng)用較多的界面捕捉方法包括相場法[19]、流體體積法[20]、擴散界面法[21]、格子-玻爾茲曼方法[22]以及水平集法[23]等. 鄭志偉等[24]基于耦合水平集法和流體體積法的CLSVOF方法建立了中空液滴撞擊固體壁面的數(shù)值模型, 研究了撞擊速度和壁面浸潤性對中空液滴撞擊壁面動力學(xué)行為和傳熱特性的影響. 李玉杰等[25]基于格子玻爾茲曼方法模擬了液滴撞擊圓柱內(nèi)壁面的過程, 分析了液滴的物性和圓柱內(nèi)壁面浸潤性等因素對撞擊行為的具體影響. 高亞軍等[26]采用水平集方法模擬了雙液滴同時撞擊固體壁面的行為, 分析了不同撞擊參數(shù)下中心射流高度和水平鋪展半徑隨時間的演化規(guī)律. 在非牛頓流體液滴方面, 基于擴散界面方法, 韓丁丁等[27]研究了用Oldroyd-B模型描述的黏彈性流體液滴撞擊疏水壁面的動力學(xué)行為,他們的數(shù)值模擬結(jié)果指出流體黏彈性對液滴的鋪展過程沒有顯著影響, 界面及移動接觸線處的彈性應(yīng)力是液滴回彈行為受到抑制的主要原因. 通過VOF方法捕捉相界面的移動, 采用Herschel-Bulkley模型耦合流體的屈服應(yīng)力-剪切變稀特性, Kim和Baek[28]研究了屈服應(yīng)力-剪切變稀液滴撞擊固體壁面的行為, 他們的研究結(jié)果表明冪律指數(shù)、非牛頓雷諾數(shù)以及韋伯數(shù)決定著撞擊液滴的鋪展階段,而在液滴的回縮階段, 非牛頓毛細數(shù)起著主要作用. 通過相場法捕捉相界面的移動, 采用Giesekus模型耦合流體的流變學(xué)特性, Wang等[29]研究了黏彈性-剪切變稀液滴撞擊固體壁面的行為, 他們指出慣性力、流體的黏性以及彈性共同影響著液滴撞擊壁面后的運動狀態(tài).

上述的分析表明, 目前關(guān)于剪切變稀流體液滴撞擊固體壁面行為的有限研究集中于實驗工作, 尚無學(xué)者對這一在相關(guān)領(lǐng)域中有極高應(yīng)用價值的行為開展系統(tǒng)的數(shù)值模擬研究. 本文基于有限元耦合水平集方法捕捉相界面的移動, 通過修正的冪律模型描述流體的剪切變稀特性, 構(gòu)建了剪切變稀液滴撞擊固體壁面的數(shù)值模型, 深入分析了剪切變稀流體液滴撞擊固體壁面過程中各參數(shù)的影響機制.

2 物理數(shù)學(xué)模型

2.1 控制方程

目前采用水平集方法求解界面演化尤其是涉及拓撲結(jié)構(gòu)變化的界面演化問題已發(fā)展的較為成熟, 本研究采用水平集方法捕捉氣液兩相界面的移動, 用水平集方程中的函數(shù)值來表示出具有限定厚度的液滴界面. 從相界面一側(cè)過渡到另一側(cè)時,值從0到1逐漸變化, 本文中純氣相中, 液相中相界面的位置由的等值線表示.通過求解水平集方程中的值來獲取兩相界面的位置. 本文求解的水平集方程為

采用界面捕捉法求解相界面移動的困難在于使高密度與黏度比的兩相在界面處平滑過渡, 本文采用如下的Heaviside光滑函數(shù)對相界面進行光滑處理,

光滑處理后, 界面處流體的密度和黏度可表示為

Navier-Stokes (N-S)方程可以用于描述不可壓縮流體質(zhì)量和動量的傳輸特性. 在固定的歐拉坐標系中, 連續(xù)性方程以及不可壓縮, 耦合界面張力和重力作用的N-S方程如下:

為了耦合剪切變稀特性對液滴撞擊固體壁面后的影響, 數(shù)值模型中采用冪律模型描述液滴剪切黏度隨剪切速率的動態(tài)變化,

式中,k是稠度指數(shù),m是冪律指數(shù). 同時為避免冪律模型在剪切速率過低時黏度失真的情況, 將冪律模型通過下式[30]進行修正:

式中D是變形張量.

通過添加壁面作用力的形式[31]耦合了壁面浸潤性對液滴鋪展過程的影響, 作用于壁面處流體的力

2.2 數(shù)值模型驗證

本文建立二維軸對稱模型模擬液滴撞擊壁面的動態(tài)過程. 常見的初始液滴位置設(shè)置方式有兩種: 初始液滴與壁面相切或初始液滴置于壁面上方并與壁面保持一段距離, 本文采取第二種設(shè)置形式, 采用此種設(shè)置方式考慮了液滴撞擊壁面過程中液滴下方氣體的影響, 更為貼近真實情況. 圖1為數(shù)值模型的計算域示意圖, 液滴關(guān)于圖左側(cè)的對稱軸對稱, 在垂直壁面的初始速度的作用下撞擊壁面并在壁面上鋪展. 其中, 上邊界及右邊界均為法向應(yīng)力等于零的開邊界, 下邊界為浸潤壁邊界條件. 對所有計算區(qū)域采用如圖2所示的自由四邊形網(wǎng)格進行剖分.

圖1 數(shù)值模擬計算域設(shè)置Fig. 1. Computation domain for numerical simulation.

圖2 計算域網(wǎng)格剖分Fig. 2. Grid generation of computation domain.

為了系統(tǒng)地分析各因素對液滴撞擊固體壁面后鋪展行為的影響, 本文引入如下的無量綱參數(shù)討論液滴的撞擊鋪展結(jié)果, 包括韋伯數(shù)非牛頓雷諾數(shù)無量綱時間無量綱直徑無量綱高度其中,為液滴在豎直方向上的速度分量,為液滴初始直徑,為液滴鋪展過程中的動態(tài)直徑,為液滴鋪展過程中的動態(tài)高度,t為時間.

表1 數(shù)值模擬參數(shù)設(shè)置Table 1. Symbols and constants in numerical simulation.

圖3 本文數(shù)值模擬與Lim等[33]相場模擬的結(jié)果比較Fig. 3. Comparisons between simulation results in this work and phase-field simulation results[33].

為了驗證數(shù)值模型的準確性, 基于建立的數(shù)值模型, 首先模擬了純水液滴撞擊固體壁面的情況,并將數(shù)值計算結(jié)果與Lim等[33]的數(shù)值計算結(jié)果進行了對比驗證. 模擬的相關(guān)參數(shù)設(shè)置見表1, 分析圖3中液滴撞擊固體壁面后無量綱直徑隨無量綱時間的變化, 可以發(fā)現(xiàn)本文的模擬結(jié)果與Lim等[33]的模擬結(jié)果的誤差較小(誤差這說明了本文建立的數(shù)值模型能夠較準確地模擬液滴撞擊固體壁面的行為. 同時為了驗證數(shù)值模型預(yù)測非牛頓流體撞擊壁面行為的準確性, 基于建立的模型, 模擬了稠度系數(shù)為 0.208 Pa·sn, 冪律指數(shù)為 0.4的剪切變稀液滴撞擊疏水壁面的行為, 并將模擬結(jié)果與實驗結(jié)果[11]進行了對比. 圖4為文獻[11]中實驗結(jié)果與本文模型模擬結(jié)果的對比, 可以看出不同時刻的數(shù)值計算結(jié)果均與實驗結(jié)果高度吻合, 說明了本文建立的數(shù)值模型預(yù)測非牛頓流體液滴撞擊壁面行為的準確性. 兩相流模擬對網(wǎng)格的精度要求較高, 網(wǎng)格的尺寸大小直接影響著數(shù)值結(jié)果的精度,圖5為本研究的網(wǎng)格尺寸無關(guān)性驗證, 所選擇的模擬參數(shù)與表1相同. 考察3種不同的網(wǎng)格數(shù)量, 依次為 23221, 51203, 90724. 選擇直線 (坐標 (0, 0)至(0, 10))上無量綱速度的徑向分量作為考察對象, 從圖5可以看出網(wǎng)格數(shù)量51203和90724的數(shù)值結(jié)果基本重合, 網(wǎng)格收斂性較好, 因此以下研究的網(wǎng)格數(shù)量均為51203.

圖4 本文數(shù)值模擬與 German 等[11]的實驗結(jié)果比較 (左側(cè)為實驗結(jié)果, 右側(cè)為模擬結(jié)果)Fig. 4. Comparisons between numerical simulation in this work and experiment results[11], the left half of each image is obtained from the experiment, while the right half is the snapshots from our simulation.

圖5 液滴撞擊壁面的網(wǎng)格收斂性驗證Fig. 5. Mesh convergence study of droplet impact on solid surfaces.

圖6 不同冪律指數(shù)流體剪切黏度隨剪切速率的變化Fig. 6. Variation of shear viscosity with shear rate at different power-law index.

3 結(jié)果分析與討論

3.1 剪切變稀特性的影響

為了研究剪切變稀特性對液滴撞擊固體壁面后鋪展行為的影響, 在不改變表1中除黏度以外參數(shù)的情況下, 對比研究了牛頓流體液滴和三種表現(xiàn)出不同剪切變稀程度液滴的撞擊壁面行為. 圖6為模擬研究中不同冪律指數(shù)流體剪切黏度隨剪切速率的變化, 稠度系數(shù)和冪律指數(shù)的具體數(shù)值選擇參考了Lindner等[34]對黃原膠稀溶液的流變學(xué)研究.由冪律模型的本構(gòu)關(guān)系可知在低剪切速率時, 其剪切黏度會趨近一個無窮大的值, 所以本文在剪切速率為時將曲線截斷, 對黏度進行修正.

圖7 不同冪律指數(shù)時無量綱直徑隨無量綱時間的變化(We = 4.53)Fig. 7. Dimensionless diameter of droplet spread varying with dimensionless time at different power-law index (We =4.53).

圖8 不同冪律指數(shù)時無量綱高度隨無量綱時間的變化(We = 4.53)Fig. 8. Dimensionless height of droplet spread varying with dimensionless time at different power-law index (We =4.53).

從圖7和圖8可以發(fā)現(xiàn), 對于表現(xiàn)出剪切變稀特性的液滴, 當m較小時它的鋪展速度更快(更大的斜率), 但由于其在壁面上的鋪展直徑更大, 所以達到最大鋪展的時間仍隨著m減小而增大. 雖然鋪展直徑的增大會導(dǎo)致液滴與壁面間的摩擦能量耗散增加, 但從鋪展直徑的變化曲線可以看出, 黏性耗散在此過程中起著更大的影響. 液滴達到最大鋪展直徑后的回縮階段已經(jīng)不受慣性力的影響, 此時主導(dǎo)液滴運動的主要是表面張力和液滴黏性力,同樣由于剪切變稀特性的存在, 液滴的回縮速度更快 (更大斜率絕對值).

圖9對比了m= 0.85和0.80的液滴撞擊壁面(接觸角為55°)時液滴形貌與速度矢量隨時間的變化情況.時, 可以發(fā)現(xiàn)m= 0.80的液滴出現(xiàn)內(nèi)部薄液膜斷裂的情況, 液滴內(nèi)部液膜斷裂的同時伴隨著一個中心小液滴的生成. 在液滴的回縮過程中, 中心小液滴會重新合并到主液滴中. 這種液滴鋪展過程中出現(xiàn)內(nèi)部液膜斷裂的行為與前人[35,36]在數(shù)值模擬研究中發(fā)現(xiàn)的行為類似, 本文中出現(xiàn)此現(xiàn)象的原因是m= 0.80液滴的表觀黏度最小, 液滴撞擊壁面過程中的黏性耗散較小, 從圖7中液滴內(nèi)部黏度場分布也可以看出, 相較于冪律指數(shù)較大的液滴, 此時液滴內(nèi)部液膜處的剪切較強,剪切黏度較小, 抵抗形變的能力最弱, 最終在表面張力和慣性力的驅(qū)使下, 內(nèi)部液膜斷裂.

圖9 不同冪律指數(shù)時液滴撞擊壁面過程(a) m = 0.85, Ren = 24.37, We = 4.53; (b) m = 0.80, Ren = 29.50, We = 4.53Fig. 9. Process of droplet impact on surface at different power-law index: (a) m = 0.85, Ren = 24.37, We = 4.53; (b) m = 0.80, Ren =29.50, We = 4.53.

3.2 壁面浸潤性的影響

液滴撞擊壁面后的鋪展行為被氣液固三相耦合作用, 固體壁面的浸潤性直接影響著壁面與液滴間力的相互作用. 通常采用接觸角衡量壁面的浸潤程度, 3.1節(jié)中已經(jīng)研究了接觸角為55°情況下液滴撞擊固體壁面時的動力學(xué)行為, 為了研究壁面浸潤性對液滴撞擊壁面后鋪展行為的影響, 繼續(xù)模擬了壁面接觸角分別設(shè)置為100°和160°時液滴撞擊固體壁面的行為, 探討壁面浸潤性對牛頓流體液滴和剪切變稀液滴撞擊壁面行為的影響差異, 圖10和圖11為數(shù)值計算結(jié)果的分析圖.

從圖10和圖11可以看出, 隨著壁面浸潤性的減小 (接觸角增大), 液滴撞擊壁面后的最大無量綱直徑增大, 但壁面浸潤性對剪切變稀液滴的最大無量綱直徑影響更為顯著. 接觸角為55°情況下,牛頓流體液滴撞擊壁面時, 液滴在達到最大鋪展后基本不表現(xiàn)出回縮行為; 接觸角為100°時液滴撞擊壁面后鋪展范圍減小, 與壁面摩擦導(dǎo)致的能量耗散較小, 達到最大鋪展后, 液滴的表面張力和壁面力驅(qū)使液滴回縮, 但此時并未彈起.接觸角為160°時, 液滴撞擊過程中的摩擦耗散能量繼續(xù)減少, 由(11)式可知壁面對液滴的作用力繼續(xù)增大,液滴最終在回縮過程中彈起.

圖10 不同冪律指數(shù)時無量綱直徑隨無量綱時間的變化(We = 4.53, Ren = 13.75 (m = 1.00), Ren = 24.37 (m =0.85))Fig. 10. Dimensionless diameter of droplet spread varying with dimensionless time at different power-law index (We =4.53, Ren = 13.75 (m = 1.00), Ren = 24.37 (m = 0.85)).

圖11 不同冪律指數(shù)時無量綱高度隨無量綱時間的變化(We = 4.53, Ren = 13.75 (m = 1.00), Ren = 24.37 (m =0.85))Fig. 11. Dimensionless height of droplet spread varying with dimensionless time at different power-law index (We =4.53, Ren = 13.75 (m = 1.00), Ren = 24.37 (m = 0.85)).

圖12 (a)給出了接觸角為160°時牛頓流體液滴撞擊壁面過程中內(nèi)部速度矢量隨無量綱時間的變化. 在液滴的撞擊瞬間, 液滴上部的速度幾乎與撞擊前一致, 此時在液滴近壁中心處的速度接近于零, 在液滴近壁邊緣處則會產(chǎn)生較大的橫向速度.隨著液滴在壁面上的鋪展, 液滴邊緣的橫向速度逐漸減小, 液滴達到最大鋪展時, 液滴邊緣的橫向速度減小為零. 在表面張力的作用下, 液滴開始回縮,此時液滴近壁中心處的速度依然接近于零, 液滴上部的速度矢量方向向下, 但液滴邊緣出現(xiàn)指向液滴中心的速度 (t= 0.64), 這導(dǎo)致此時液滴內(nèi)部出現(xiàn)了較小的速度環(huán)流. 隨著液滴的回縮, 液滴橫向回縮速度逐漸增大, 在液滴即將脫離壁面回彈起的時刻, 近壁中心處的橫向速度達到最大值, 液滴脫離壁面后內(nèi)部基本不存在徑向的速度分量. 圖12(b)為剪切變稀液滴撞擊壁面過程中內(nèi)部速度矢量分布的變化情況, 剪切變稀液滴撞擊壁面后其內(nèi)部的速度矢量分布變化與牛頓流體基本一致, 但由于流體的剪切變稀特性, 液滴鋪展過程中的黏性耗散小于牛頓流體液滴, 黏性耗散的減少使得更多的慣性能量轉(zhuǎn)化為表面能量, 在表面張力的驅(qū)使下液滴回縮的速度更快, 液滴彈起壁面的時間略早于牛頓流體液滴, 最終回彈高度也略高.

圖12 不同冪律指數(shù)時液滴撞擊壁面過程(a) m = 1, Ren = 13.75, We = 4.53; (b) m = 0.80, Ren = 29.50, We = 4.53Fig. 12. Process of droplet impact on surface at different power-law index: (a) m = 1, Ren = 13.75, We = 4.53; (b) m = 0.80, Ren =29.50, We = 4.53.

3.3 最大無量綱直徑預(yù)測模型

由圖13可知(12)式的模型預(yù)測結(jié)果與本文數(shù)值計算結(jié)果具有較高的一致性.

表2 最大無量綱直徑預(yù)測模型Table 2. Prediction models of maximum dimensionless factor.

圖13 最大無量綱直徑模型預(yù)測值與數(shù)值計算值對比Fig. 13. Comparison of between model prediction values and simulation data.

4 結(jié)論

本文基于有限元法, 通過水平集方法捕捉相界面的移動, 構(gòu)建了液滴撞擊固體壁面的數(shù)值模型.采用修正的冪律模型描述流體的剪切變稀特性, 通過研究液滴形貌及無量綱參數(shù)隨無量綱時間的變化規(guī)律, 分析了剪切變稀特性對液滴撞擊不同浸潤性壁面行為的影響機制, 得到如下結(jié)果.

1)液滴撞擊壁面的過程中, 剪切變稀液滴內(nèi)部剪切速率的改變使液滴的黏度動態(tài)變化. 冪律指數(shù)m較小的液滴, 液滴形貌變化更為顯著. 在接觸角為55°時液滴撞擊壁面的情況下, 當冪律指數(shù)m小于0.85時, 液滴撞擊壁面后開始表現(xiàn)出振蕩行為, 當冪律指數(shù)m小于0.80時, 液滴撞擊壁面過程中出現(xiàn)內(nèi)部液膜斷裂的情況.

2)壁面浸潤性對剪切變稀液滴撞擊壁面行為的影響更為顯著, 在接觸角為100°時液滴撞擊壁面情況下,m= 0.85的剪切變稀流體表現(xiàn)出區(qū)別與牛頓流體的振蕩行為. 當接觸角為160°時, 牛頓流體液滴和剪切變稀液滴撞擊壁面時均會彈起壁面,但因黏性耗散較小, 剪切變稀液滴彈起的速度更快.

3)基于數(shù)值計算結(jié)果和前人的經(jīng)典黏性耗散模型, 本文提出了接觸角為55°時剪切變稀液滴撞擊壁面后最大無量綱直徑的預(yù)測模型, 模型預(yù)測結(jié)果與本文的數(shù)值計算結(jié)果取得了較高的一致性.