剪切變稀液滴撞擊不同浸潤性壁面的數值模擬研究*

沈學峰 曹宇 王軍鋒 劉海龍

(江蘇大學能源與動力工程學院, 鎮江 212013)

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

1 引言

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

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

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

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

2 物理數學模型

2.1 控制方程

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

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

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

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

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

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

式中D是變形張量.

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

2.2 數值模型驗證

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

圖1 數值模擬計算域設置Fig. 1. Computation domain for numerical simulation.

圖2 計算域網格剖分Fig. 2. Grid generation of computation domain.

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

表1 數值模擬參數設置Table 1. Symbols and constants in numerical simulation.

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

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

圖4 本文數值模擬與 German 等[11]的實驗結果比較 (左側為實驗結果, 右側為模擬結果)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 液滴撞擊壁面的網格收斂性驗證Fig. 5. Mesh convergence study of droplet impact on solid surfaces.

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

3 結果分析與討論

3.1 剪切變稀特性的影響

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

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

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

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

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

圖9 不同冪律指數時液滴撞擊壁面過程(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節中已經研究了接觸角為55°情況下液滴撞擊固體壁面時的動力學行為, 為了研究壁面浸潤性對液滴撞擊壁面后鋪展行為的影響, 繼續模擬了壁面接觸角分別設置為100°和160°時液滴撞擊固體壁面的行為, 探討壁面浸潤性對牛頓流體液滴和剪切變稀液滴撞擊壁面行為的影響差異, 圖10和圖11為數值計算結果的分析圖.

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

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

圖12 不同冪律指數時液滴撞擊壁面過程(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 最大無量綱直徑預測模型

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

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

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

4 結論

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

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

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

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