方可寧，單彥廣，袁俊杰

（上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093）

0 引 言

固著液滴蒸發(fā)是指附著在固體壁面上的液滴蒸發(fā)過程。隨著蒸發(fā)進(jìn)行，液滴體積不斷減小，同時(shí)液滴外部形狀不斷發(fā)生變化，不飽和蒸汽發(fā)生擴(kuò)散。液滴蒸發(fā)過程涉及復(fù)雜的傳熱傳質(zhì)，受周圍環(huán)境、基板溫度、壁面潤濕性以及液滴特性等因素影響。許多學(xué)者通過實(shí)驗(yàn)、理論以及數(shù)值模擬等方法做了大量研究工作。PICKNETT等[1]對(duì)固著在基板上的甲基乙酰乙酸液滴緩慢蒸發(fā)過程的質(zhì)量及體積變化進(jìn)行了研究，得出固著液滴蒸發(fā)過程中的兩種模式：接觸角不變模式和接觸半徑不變模式。BIRDI等[2]對(duì)玻璃基板上不同尺寸的純水液滴和辛烷液滴的蒸發(fā)過程中液滴質(zhì)量和接觸角的變化情況進(jìn)行研究，結(jié)果表明液滴蒸發(fā)速率由固著液滴的初始接觸半徑?jīng)Q定。此外，許多學(xué)者對(duì)親水[3-5]或疏水[6-8]基板上的液滴蒸發(fā)過程中接觸角、接觸半徑以及形狀變化也進(jìn)行了研究。由于實(shí)驗(yàn)研究難以通過對(duì)蒸發(fā)通量及蒸發(fā)速率的精確控制來總結(jié)液滴的傳熱規(guī)律，因此研究人員采用數(shù)值模擬的方法對(duì)蒸發(fā)相變過程進(jìn)行研究。王寶和等[9]研究發(fā)現(xiàn)，隨液滴蒸發(fā)過程的進(jìn)行，因基板與液滴間的熱量傳遞、蒸汽向大氣環(huán)境中輸送、液滴質(zhì)量減少等因素，引起氣液界面的變化，使得液滴蒸發(fā)過程內(nèi)部流動(dòng)復(fù)雜[10-11]。DUH等[12]提出液滴緩慢蒸發(fā)的數(shù)值模型，用于研究固著在等溫基板上的液滴在蒸發(fā)過程中內(nèi)部流場(chǎng)的分布情況。文章指出液滴內(nèi)部的運(yùn)動(dòng)是由于流場(chǎng)內(nèi)存在浮力和毛細(xì)對(duì)流，從而導(dǎo)致流場(chǎng)呈現(xiàn)多種結(jié)構(gòu)。在此基礎(chǔ)上，MARZO等[13]提出將一個(gè)半無限大固體與液滴耦合，利用有限體積法進(jìn)行數(shù)值模擬，其中假定液滴在蒸發(fā)過程中保持球冠型，且液滴的體積變化與蒸發(fā)時(shí)間呈指數(shù)關(guān)系。在此基礎(chǔ)上，饒玲等[14]將橢球冠模型用于模擬液滴蒸發(fā)過程，比對(duì)實(shí)際實(shí)驗(yàn)結(jié)果[1]，相對(duì)于前人使用的球冠模型，橢球冠模型能更準(zhǔn)確地描述實(shí)際的物理過程。

相對(duì)于目前普遍應(yīng)用的計(jì)算流體力學(xué)方法（computational fluid dynamics, CFD），格子玻爾茲曼方法（Lattice Boltzmann Method, LBM）作為介觀方法，由于其邊界處理簡便，不需要額外追蹤界面等優(yōu)點(diǎn)，已應(yīng)用于不少相變傳熱傳質(zhì)問題研究。已經(jīng)有學(xué)者建立了多種LBM模型，包括顏色流體模型[15-16]、偽勢(shì)模型[17-18]、自由能模型[19]和動(dòng)力學(xué)模型[20-21]。LEE等[22]提出不可壓縮二元流體的邊界條件改進(jìn)，將模擬液滴撞擊干燥表面的過程與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，結(jié)果表明對(duì)于中等接觸角的液滴能準(zhǔn)確預(yù)測(cè)，但對(duì)大接觸角的液滴則沒有那么準(zhǔn)確。此后，LATIFIYAN等[23]將二維流體模型與熱格子結(jié)合，模擬了多孔介質(zhì)中液體滲透過程，研究無量綱參數(shù)、孔隙度、接觸角和密度比對(duì)穿透率及蒸發(fā)速率的影響。GONG等[24]對(duì)熱格子模型進(jìn)行改進(jìn)，改進(jìn)的模型可有效模擬多相流動(dòng)和相變情況，而且能解決更大溫度范圍內(nèi)的復(fù)雜流體問題。LI等[25]基于現(xiàn)有的蒸發(fā)相變模型和有焓變的固液相變模型，建立新模型模擬干燥飽和蒸汽瞬時(shí)冷凝后在傾斜的疏水基板上的凍結(jié)過程，指出重力和基板潤濕性對(duì)液滴形態(tài)有重要的影響。

針對(duì)目前基于CFD方法模擬液滴蒸發(fā)大多忽略重力的影響，將液滴蒸發(fā)過程的外形近似為球冠形的局限，采用改進(jìn)的熱格子玻爾茲曼方法，研究重力場(chǎng)下基板潤濕性和初始環(huán)境溫度對(duì)純液滴蒸發(fā)的影響，并考慮液滴在滴落基板后因與基板之間親疏水性不同而導(dǎo)致的鋪展或收縮過程，預(yù)測(cè)液滴蒸發(fā)過程內(nèi)部流場(chǎng)、外形演變、接觸直徑以及剩余質(zhì)量的變化。

1 數(shù)值模擬

1.1 多相流LBM模型

LBM通過流體粒子分布函數(shù)描述流場(chǎng)內(nèi)固定格子的運(yùn)動(dòng)過程，演化方程為：

其中：x是計(jì)算流場(chǎng)內(nèi)的一個(gè)格點(diǎn)；t為當(dāng)前時(shí)間步；ei是流體粒子的離散速度；δt是離散的時(shí)間步長；τf為無量綱松弛時(shí)間；fi為流體粒子的速度分布函數(shù)。此外，為平衡態(tài)分布函數(shù)，其表達(dá)式為：

式中：u為格子速度；ωi為加權(quán)系數(shù)；cs為格子聲速，在二維D2Q9模型中取值離散速度ei具體公式如下：

式中：Δu是在δt時(shí)刻由于力F導(dǎo)致的速度差，表示為其中，力項(xiàng)F表示為粒子點(diǎn)受到各項(xiàng)力的總和：

Fint（x）為各粒子間的相互作用力，表示為：

ψ（x）為粒子間的相互作用勢(shì)，主要由粒子的密度ρ及壓力p決定：

其中，c0為作用勢(shì)系數(shù)，取值c0=6。

此外，β作為權(quán)重因子，用于提高數(shù)值計(jì)算的穩(wěn)定性，根據(jù)狀態(tài)方程取值不同，取值β=1.16。在二維D2Q9模型中，

Fs（x）為流體與固壁面之間的相互作用力，表達(dá)式如下：

其中：gs表示流固相互作用強(qiáng)度的參數(shù)，用于控制基板的潤濕條件；s（x）用于區(qū)別固-流兩相，流體相s（x）=0，固體相s（x）=1。

Fg（x）為重力，可表示為：

其中：g為重力加速度。

實(shí)際粒子的宏觀密度及速度由下式得到：

但式（12）得到的速度并不是真正的實(shí)際速度，實(shí)際速度U需要考慮力F（x）的作用，表達(dá)式為：

由于模擬兩相間的密度比較大，因此選用R-K狀態(tài)方程能較好地模擬物理過程，其方程如下：

1.2 能量方程模型

溫度分布函數(shù)描述流場(chǎng)內(nèi)溫度場(chǎng)的變化過程，演化方程為：

其中：τT為無量綱松弛時(shí)間；gi為流體溫度分布函數(shù)；為參與氣液相變的源項(xiàng)。

溫度的平衡態(tài)分布函數(shù)[25]：

其中，相變?cè)错?xiàng)由文獻(xiàn)[24]得出，

溫度可由統(tǒng)計(jì)計(jì)算得到：

1.3 模型驗(yàn)證

表1和表2給出了數(shù)值模擬中的格子單位和實(shí)際對(duì)應(yīng)的參數(shù)值。如無特別說明，下文變量單位為表1的格子單位。

表1 格子玻爾茲曼方法中單位[26]Table 1 Units in Lattice Boltzmann Method

表2 模擬計(jì)算中參數(shù)值Table 2 Parameter value in simulation

1.3.1 Laplace 驗(yàn)證

模擬計(jì)算域在Lx×Ly內(nèi)半徑為R的懸浮液滴在環(huán)境溫度恒定為Ts=0.9×Tc條件下最終穩(wěn)定的問題。根據(jù)Maxwell重構(gòu)方程可得，在Ts=0.9×Tc的環(huán)境溫度下，液相初始密度為ρl=5.9079，氣相密度為ρg=0.5801。為保證模擬結(jié)果有效性，分別取液滴初始半徑為24、28、32、36和40的情況，進(jìn)行模擬。計(jì)算域四周的邊界條件均為周期邊界條件。懸浮液滴在表面張力σ作用下保持球冠形，達(dá)到穩(wěn)定后，求得液滴界面內(nèi)外的壓力躍變差值Δp=pl-pv的數(shù)值與半徑的倒數(shù)1/R的關(guān)系。結(jié)果如圖1所示，Δp=pl-pv與1/R呈線性關(guān)系，符合Laplace定律。

圖1 Laplace驗(yàn)證Fig.1 Laplace validation

1.3.2 網(wǎng)格無關(guān)性驗(yàn)證

模擬計(jì)算域?yàn)長x×Ly內(nèi)固著液滴在加熱基板上受重力作用的蒸發(fā)過程。同樣根據(jù)Maxwell重構(gòu)方程分別求出不同飽和溫度下的氣液兩相臨界密度，即對(duì)應(yīng)初始密度的取值，基板潤濕性為親水，接觸角θ=60°，即三相接觸線所在的氣-液界面的切線與基板間的夾角為60°，通過流固相互作用強(qiáng)度參數(shù)與接觸角的擬合，即取重力加速度取值選用不同大小的網(wǎng)格數(shù)，對(duì)液滴蒸發(fā)過程中基板上液滴的高H與接觸直徑D之間比值的變化進(jìn)行記錄。每10 000次迭代計(jì)算輸出一次結(jié)果，數(shù)據(jù)匯總見表3。

表3 網(wǎng)格無關(guān)性驗(yàn)證Table 3 Grid independence validation

由表3可知，不同網(wǎng)格數(shù)下的H/D值變化趨勢(shì)基本一致，隨蒸發(fā)進(jìn)行而不斷增大。網(wǎng)格數(shù)在50 × 50與100 × 100之間的H/D的最大相對(duì)差值為0.96%，網(wǎng)格數(shù)在100 × 100與200 × 200之間的相對(duì)差值最大為0.06%，而網(wǎng)格數(shù)200 × 200與400 × 400之間的最大相對(duì)誤差為0.02%。因此計(jì)算域Lx×Ly選用200 × 200網(wǎng)格數(shù)進(jìn)行模擬計(jì)算。

2 結(jié)果及分析

模擬計(jì)算域?yàn)長x×Ly，假定液滴初始滴落到基板上呈半圓形（即接觸角為90°）。受基板潤濕性影響，液滴與基板接觸后發(fā)生鋪展或收縮現(xiàn)象，穩(wěn)定后達(dá)到設(shè)定的靜態(tài)接觸角，然后開始蒸發(fā)直至蒸發(fā)完成。為簡化計(jì)算，假定流體物性參數(shù)在蒸發(fā)過程中保持不變，忽略靜壓力對(duì)液滴蒸發(fā)的影響，加熱基板為光滑基板，具體物理模型如圖2所示。Lx×Ly恒定為200 × 200網(wǎng)格數(shù)，藍(lán)色區(qū)域?yàn)槌跏及雸A形液滴，放置在加熱基板的中心位置，初始半徑R=40，加熱基板溫度為Tw=0.9×Tc，初始環(huán)境溫度Ts=0.85Tc，液滴初始溫度與環(huán)境溫度相等。據(jù)初始環(huán)境溫度及Maxwell重構(gòu)方程，氣相及液相密度分別取ρg=0.5801、ρl=5.9079。該模擬計(jì)算域左右邊界設(shè)置為周期性邊界條件，系統(tǒng)為開放系統(tǒng)，因此上邊界為對(duì)流邊界條件，下邊界出口為加熱基板，設(shè)定下邊界溫度為加熱溫度，設(shè)置為無滑移邊界條件。

圖2 計(jì)算模型Fig.2 Calculation model

為了研究蒸發(fā)過程液滴形狀的變化以及流體的傳熱流動(dòng)特性，通過改變液滴重力及基板潤濕性、初始環(huán)境溫度等條件，模擬不同條件下的液滴蒸發(fā)過程，考察液滴接觸直徑、蒸發(fā)速率以及蒸發(fā)壽命的變化情況。其中液滴接觸直徑，即液滴與基板間的三相接觸線直徑，記為D。蒸發(fā)速率即計(jì)算域內(nèi)液相質(zhì)量變化量與初始液相總質(zhì)量之比記為Vt，其中液相的定義為密度值的區(qū)域。蒸發(fā)壽命為液滴在加熱基板上完全蒸干所需要的時(shí)間。對(duì)演化時(shí)間進(jìn)行無量綱化，特征時(shí)間定義為同樣，文中模擬變量均為LBM對(duì)應(yīng)的無量綱單位。

2.1 重力對(duì)液滴蒸發(fā)過程的影響

為考察重力對(duì)液滴在加熱基板上鋪展及蒸發(fā)過程的影響，模擬考慮重力和忽略重力兩種情況的鋪展蒸發(fā)過程?；鍧櫇裥詾橛H水，接觸角θ=60°，流固相互作用強(qiáng)度參數(shù)和重力加速度取值同1.3.2節(jié)。

圖3為液滴蒸發(fā)過程演化圖。由圖可知，兩種情況下，受基板潤濕性影響，液滴放置在基板直至穩(wěn)定階段，三相接觸線均會(huì)發(fā)生鋪展現(xiàn)象，但重力情況下的鋪展現(xiàn)象更明顯。穩(wěn)定后液滴進(jìn)行蒸發(fā)，液滴體積隨蒸發(fā)時(shí)間不斷減小，同時(shí)接觸角基本不變，直至蒸發(fā)完成。相同時(shí)間內(nèi)，重力情況下液滴蒸發(fā)更快，這是由于液滴受重力作用，液滴與加熱基板的接觸面積增大，通過接觸面?zhèn)鬟f至液滴內(nèi)部的熱量也就更多，因此液滴在重力條件下能較快蒸干。

圖3 液滴蒸發(fā)過程演化圖Fig.3 Evolution diagram of droplet evaporation process

為驗(yàn)證上述結(jié)論，圖4給出了液滴鋪展及蒸發(fā)過程中接觸直徑隨時(shí)間的變化趨勢(shì)。由圖4可知，液滴放置在基板后，受基板潤濕性影響，接觸直徑逐漸增大至穩(wěn)定。隨后蒸發(fā)進(jìn)行，液滴三相接觸線不斷收縮，接觸直徑隨之減小，接觸直徑與時(shí)間呈冪指數(shù)減小關(guān)系，該現(xiàn)象與文獻(xiàn)[27]結(jié)果一致。重力條件下，液滴鋪展現(xiàn)象更明顯，對(duì)應(yīng)穩(wěn)定后的接觸直徑較大。與文獻(xiàn)[27]結(jié)果不同的是，重力情況下液滴接觸直徑減小速度較快，蒸發(fā)后期同一時(shí)刻液滴的接觸直徑相對(duì)較小，完成蒸發(fā)所用的時(shí)間也較短。

圖4 液滴接觸直徑隨時(shí)間變化趨勢(shì)圖Fig.4 Trend chart of contact diameter change in droplet evaporation process

圖5給出了蒸發(fā)過程中液滴剩余質(zhì)量隨時(shí)間的變化趨勢(shì)。由圖可知，兩種重力條件下，液滴剩余質(zhì)量變化趨勢(shì)基本一致，均隨蒸發(fā)時(shí)間呈線性減小趨勢(shì)。相對(duì)于無重力條件下的剩余質(zhì)量變化，有重力情況下的斜率更大，即蒸發(fā)量較大，因而蒸發(fā)所需時(shí)間較短，能夠較早完成蒸發(fā)。

圖5 重力對(duì)液滴剩余質(zhì)量隨時(shí)間變化的影響Fig.5 Influence of gravity on the residual mass of droplets with time

由上述結(jié)論可知，重力場(chǎng)對(duì)液滴蒸發(fā)有不可忽略的影響，也即相較于常用的球冠模型，橢球冠液滴模型與實(shí)際過程更為接近[15]。因此研究均考慮重力場(chǎng)對(duì)液滴的初始鋪展或收縮過程及蒸發(fā)的影響。

2.2 潤濕性對(duì)液滴蒸發(fā)過程的影響

基板潤濕性是影響液滴鋪展及蒸發(fā)的重要因素，本節(jié)模擬不同基板潤濕性條件下液滴鋪展蒸發(fā)過程?？紤]靜態(tài)接觸角為60°、90° 和120°，即gs取值分別為0.22、0.70、1.14。其余參數(shù)與上節(jié)相同。

圖6為不同潤濕條件下液滴鋪展及蒸發(fā)過程的演化。可知基板潤濕性表現(xiàn)為親水時(shí)，液滴放置在基板后會(huì)發(fā)生鋪展現(xiàn)象，直到接觸角接近60° 時(shí)達(dá)到穩(wěn)定。液滴蒸發(fā)過程，鋪展穩(wěn)定后的接觸角保持不變，接觸直徑不斷減小，直至液滴完全蒸干。當(dāng)接觸角為90° 時(shí)，液滴維持初始形狀，不發(fā)生鋪展。蒸發(fā)過程中液滴接觸角不變，接觸直徑不斷減小，直至液滴完全蒸干。當(dāng)基板潤濕性表現(xiàn)為疏水時(shí)，液滴在滴至基板后會(huì)出現(xiàn)收縮現(xiàn)象，直至液滴形狀維持在接觸角120° 時(shí)達(dá)到穩(wěn)定。隨后液滴繼續(xù)蒸發(fā)，直至完全蒸干。比較三種基板接觸角條件下的蒸發(fā)過程可以發(fā)現(xiàn)，相同時(shí)刻下，隨著基板疏水性增強(qiáng)，殘余液滴質(zhì)量相對(duì)較大，即相應(yīng)液滴蒸發(fā)較慢。

圖6 不同潤濕性條件下液滴蒸發(fā)過程演化圖Fig.6 Evolution diagram of droplet evaporation process under different wettability conditions

為驗(yàn)證以上結(jié)論，考察基板接觸角為60°、75°、90°、105° 和120° 五種情況下液滴的蒸發(fā)壽命，如圖7。可知加熱基板接觸角越小，即越親水，液滴蒸發(fā)壽命越短；基板接觸角越大，液滴蒸發(fā)壽命相對(duì)較長?；鍨槭杷畷r(shí)，液滴的蒸發(fā)壽命隨潤濕性變化幅度更大。

為進(jìn)一步分析基板潤濕性對(duì)液滴壽命的影響，圖8給出了不同基板潤濕性條件下液滴接觸直徑隨時(shí)間的變化情況。可知，親水基板的液滴鋪展穩(wěn)定后的接觸直徑較大，隨蒸發(fā)進(jìn)行減小相對(duì)較快。接觸角為90° 時(shí)液滴維持原狀，不出現(xiàn)鋪展或收縮現(xiàn)象，接觸直徑隨時(shí)間變化基本呈冪指數(shù)相關(guān)，蒸發(fā)壽命相較親水基板的液滴變化不大。疏水基板上的液滴在收縮穩(wěn)定后，因接觸直徑較小，其減小趨勢(shì)較緩，蒸發(fā)壽命較長。

圖7 不同潤濕性條件下液滴蒸發(fā)壽命變化圖Fig.7 Diagram of droplet evaporation life under different wettability conditions

圖8 不同潤濕性條件下液滴接觸直徑隨時(shí)間變化趨勢(shì)圖Fig.8 Trend chart of contact diameter change with time under different wettability conditions

圖9為液滴剩余質(zhì)量隨時(shí)間的變化情況?？芍诓煌瑵櫇裥缘幕迳弦旱问Ｓ噘|(zhì)量變化趨勢(shì)基本一致，液滴初始質(zhì)量相等，剩余質(zhì)量均隨蒸發(fā)時(shí)間呈線性減小趨勢(shì)。

圖9 不同潤濕性條件下液滴剩余質(zhì)量隨時(shí)間變化圖Fig.9 Diagram of droplet evaporation rate with time under different wettability conditions

結(jié)合圖7和圖8可知，基板接觸角越小，液滴剩余質(zhì)量變化越大，也就越容易蒸干。

2.3 初始環(huán)境溫度對(duì)液滴蒸發(fā)過程的影響

為探究初始溫度場(chǎng)（環(huán)境溫度記Ts、基板溫度記Tw）對(duì)液滴鋪展及蒸發(fā)過程的影響，本節(jié)模擬研究不同環(huán)境溫度下液滴鋪展及蒸發(fā)過程。三個(gè)工況下的初始溫度設(shè)置如表4所示，gs始終取值0.70。

表4 不同工況下初始溫度Table 4 Initial temperature under different working conditions

圖10 不同環(huán)境溫度下液滴蒸發(fā)過程演化圖Fig.10 Evolution diagram of droplet evaporation process at different environmental temperatures

圖10給出了不同環(huán)境溫度下，液滴蒸發(fā)的演化過程圖。由圖可知，不同初始環(huán)境溫度下液滴的鋪展及蒸發(fā)過程趨勢(shì)一致。液滴在加熱基板上蒸發(fā)，整體形狀逐漸變小，而接觸角在蒸發(fā)過程中保持不變。初始環(huán)境溫度與基板溫度之間的溫度差越大時(shí)，在同一演化時(shí)間，基板上殘留的液滴質(zhì)量更小，蒸發(fā)量較大，即蒸發(fā)速率較快。

為進(jìn)一步考察初始環(huán)境溫度對(duì)液滴鋪展及蒸發(fā)過程中形狀變化的影響，統(tǒng)計(jì)了不同工況下液滴接觸直徑隨時(shí)間變化的趨勢(shì)，如圖11所示。由圖可知，各工況下液滴接觸直徑變化趨勢(shì)基本一致。液滴放置在基板直至穩(wěn)定階段，基本都存在不同程度的收縮過程，然后開始蒸發(fā)，其接觸直徑持續(xù)減小，直至完全蒸干。

圖11 不同環(huán)境溫度時(shí)接觸直徑隨蒸發(fā)時(shí)間變化趨勢(shì)圖Fig.11 Trend chart of contact diameter change with time under different initial temperature fields

3 總 結(jié)

采用LBM模擬了固著在加熱基板上純液滴的鋪展及蒸發(fā)過程，研究重力場(chǎng)、基板潤濕性和初始環(huán)境溫度對(duì)液滴蒸發(fā)過程形狀變化、液滴鋪展或收縮現(xiàn)象以及蒸發(fā)速率等傳熱規(guī)律的影響，結(jié)論如下：

（1）液滴蒸發(fā)過程中需考慮重力場(chǎng)的影響。重力作用下基板上液滴穩(wěn)定后呈橢球冠，相較無重力液滴鋪展面積更大，蒸發(fā)速率大，蒸發(fā)所需時(shí)間較短。

（2）基板潤濕性不同，液滴鋪展或收縮穩(wěn)定后與基板的接觸角不同。相同質(zhì)量條件下，液滴與基板接觸角越大，發(fā)生收縮現(xiàn)象穩(wěn)定后接觸直徑越小，因此蒸發(fā)速率越小，所需的蒸發(fā)時(shí)間更長。

（3）環(huán)境與基板溫度的溫差較大時(shí)，傳熱驅(qū)動(dòng)力更大，強(qiáng)化了熱量傳遞，液滴蒸發(fā)更快。