2理論模型及計算方法

2.1理論模型

假設一不可壓縮液滴在均勻加熱的固體表面上不斷蒸發(fā),如圖1所示,液滴初始最大厚度為,流動方向尺度為L,壁面溫度為w.由于ε=?1(圖1中的H?,L?僅為表示清楚,不代表真實比例),因此可應用潤滑理論(上標～表示有量綱量,下同).

圖1液滴在加熱表面鋪展示意圖

液滴運動過程的控制方程包括連續(xù)性方程、動量方程和能量方程:

界面處的蒸發(fā)流量與界面溫度有關,可表示為

為避免接觸線處出現(xiàn)應力奇點,在水平方向施加Navier滑移邊界條件:

壁面溫度滿足

由于固-液、固-氣界面張力目前仍難以通過實驗直接測量,假設固-液、固-氣界面張力不隨溫度改變,而氣-液界面張力與溫度滿足線性關系:

將控制方程和邊界條件采用下式進行無量綱化:

通過無量綱變換和保留數(shù)量級≥Ο(ε),可得無量綱控制方程組為

無量綱邊界條件為無量綱形式的界面張力與溫度關系為

液滴表面溫度為

采用積分法可得液滴厚度的演化方程為

(27)式右側依次代表毛細力、重力、熱毛細力和蒸發(fā)在液滴運動過程中的作用.

液滴在理想光滑表面上,在三相接觸線處,Young方程為

為描述接觸角遲滯影響,采用以下線性經(jīng)驗關聯(lián)式描述前進接觸角和后退接觸角與遲滯角之間的關系:

式中,A≥0是與粗糙度無關,與液體性質有關的常數(shù);遲滯角Δθ與表面粗糙度和化學異質性有關.

接觸線處的液滴高度為零,即h(x=xcl=xcr,t)=0,其中xcl和xcr分別代表左側和右側接觸線位置.采用下式描述接觸線移動速率與接觸角間的關系[28](下標c代表接觸線):

式中,B和m為常數(shù);θ為動態(tài)接觸角;υc>0代表接觸線向外鋪展,υc<0代表向內收縮.

計算初始條件為

式中,F(x)=0.5[1+tanh(20x)].

基于Karapetsas等提出的坐標變換法,將瞬態(tài)物理域(x,t)映射到固定計算域(x',t'),

液滴運動區(qū)域設定為–1≤x'≤1,時間導數(shù)表示如下:

(36)式用于替換控制方程(27)的相應項.

數(shù)值求解采用Freefem++14.3,該軟件是求解偏微分方程的強大工具,是一款免費的、開放源代碼的有限元計算軟件,目前已經(jīng)比較成熟地運用在各類數(shù)值模擬中,對于液滴、液膜的流體數(shù)值計算也很常見.該軟件采用Delaunay算法生成離散偏微分方程所需網(wǎng)格,具有網(wǎng)格自適應和移動網(wǎng)格生成的功能;采用稀疏矩陣存儲格式,內存需求少,計算速度快,是求解復雜區(qū)域問題的高效計算軟件.

2.2計算模型驗證

參照朱仙仙等用4μL去離子水在加熱溫度為60°的玻璃表面上的蒸發(fā)實驗數(shù)據(jù),換算得到無量綱參數(shù),各參數(shù)取值為:ε=0.1,A=0.39,B=0.007,m=3,C=–0.004,Bo=0.007,Ca=0.5,K=10,E=1×10-3,Ωlg=0.003.首先,為了驗證數(shù)值模擬結果的準確性,將接觸角和接觸線隨時間變化的實驗結果無量綱化后與模擬值進行對比,結果如圖2(a)和圖2(b)所示.可以看出,Δθ=0.27時的模擬值與實驗結果總體符合良好,即采用考慮接觸角遲滯的理論模型可更準確地反映液滴的運動過程.值得注意的是,朱仙仙等所給實驗結果并未呈現(xiàn)液滴鋪展階段,可能是由此階段變化較快、不易觀測記錄所致;模擬結果和實驗結果有偏差的原因可能是實驗時基底加熱時受熱不均勻或者基底傳熱能力不一致,也有可能是基于潤滑理論建立模型時,未考慮液滴蒸氣擴散所致.

圖2計算模型的驗證(a)接觸角的演化過程;(b)接觸線的演化過程;(c)網(wǎng)格無關性驗證

此外,模擬時進行了網(wǎng)格無關性驗證,在x方向上對比了網(wǎng)格數(shù)為300,500和800的結果,發(fā)現(xiàn)增加網(wǎng)格數(shù)使計算時長增大,而各特征參數(shù)計算結果數(shù)值差距不大,以圖2(c)的接觸線演化過程為例,網(wǎng)格數(shù)為500時可同時滿足精確性與高效性的要求.因此,計算中時間步長設為10–6,將[0,2]×[0,1]的計算域劃分為500×2個均勻網(wǎng)格.因方程(27)是一維的,液滴厚度h與y方向無關,所以為節(jié)約計算時間,在y方向僅劃分2層網(wǎng)格.