胡曉瑋，蔣鶴清,段淑娜，王曄春，趙 于，袁越錦*

(1.陜西科技大學(xué) 機(jī)電工程學(xué)院，陜西 西安 710021;2.西安交通大學(xué) 動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

0 引言

液滴蒸發(fā)現(xiàn)象在自然界中普遍存在，已經(jīng)受到廣泛的關(guān)注與研究，如噴墨打印[1]、DNA芯片的制造[2]、固體表面上的顆粒沉積[3]、向內(nèi)燃機(jī)中噴射燃油[4]、噴霧干燥技術(shù)[5]等.液滴蒸發(fā)過程包括液滴與外部環(huán)境的熱量交換、動(dòng)量交換及能量交換，以及液滴內(nèi)部的流動(dòng)與傳熱等，決定了相關(guān)工業(yè)工程及設(shè)備的效率，是實(shí)現(xiàn)強(qiáng)化傳熱傳質(zhì)過程的關(guān)鍵之一.

例如，燃油液滴的蒸發(fā)快慢及其與周圍環(huán)境氣體的混合程度直接決定了燃油火焰的燃燒效率、冷卻塔水滴的蒸發(fā)則決定了冷卻塔的換熱效率、化工萃取過程中工質(zhì)的蒸發(fā)決定了設(shè)備的傳質(zhì)效率等.影響液滴蒸發(fā)的因素主要有接觸角、三相接觸線、基板溫度和液滴內(nèi)部的微小流動(dòng)等，隨著研究的不斷深入與累積，液滴內(nèi)部的Marangoni流動(dòng)在液滴蒸發(fā)過程中發(fā)揮的作用正被逐步揭示.

液滴在蒸發(fā)過程中，由于液滴彎曲的自由表面及不同的熱通道長度，導(dǎo)致蒸發(fā)速率不均勻，溫度分布不均勻，在自由表面形成溫度梯度，從而引發(fā)表面張力梯度，形成Marangoni對流.Marangoni對流是一種與重力無關(guān)的、由表面引發(fā)的自然對流現(xiàn)象.Marangoni對流一般可由溫度梯度和濃度梯度兩種形式驅(qū)動(dòng)，其中溫度變化導(dǎo)致的Marangoni對流控制方式相對容易，可以較好地實(shí)現(xiàn)流動(dòng)的調(diào)控.深入研究Marangoni對流在蒸發(fā)過程中的內(nèi)在機(jī)理，有助于調(diào)控蒸發(fā)液滴的運(yùn)動(dòng)特性，進(jìn)而改善工藝，增加產(chǎn)能，創(chuàng)大收益.

早在1986年，羅銳[6]研究了平板上蒸發(fā)液滴內(nèi)部的微小流動(dòng)對蒸發(fā)過程的影響；Hegseth等[7]通過實(shí)驗(yàn)研究展示了液滴內(nèi)部的自然對流.隨后，越來越多的學(xué)者致力于更為細(xì)致、微觀的Marangoni對流研究.目前，Marangoni對流的研究方法主要有以下三種：實(shí)驗(yàn)觀測法、理論分析法、數(shù)值模擬法.

在界面可見的尺度條件下，可通過添加示蹤劑的方法來直接觀測、研究Marangoni對流[8，9].例如，Schwabe[10]選擇鋁箔作為示蹤劑，觀察到微重力環(huán)境下加熱的自由平液面上形成的Marangoni對流，研究發(fā)現(xiàn)Marangoni對流與馬蘭戈尼數(shù)Ma數(shù)密切相關(guān).但由于添加了示蹤流場的物質(zhì)，會(huì)不同程度地影響流體界面的性質(zhì)及相應(yīng)的流場分布特性.

而光學(xué)測量法具有不受示蹤劑影響、對流場介質(zhì)無干擾的特點(diǎn)，因此被廣泛應(yīng)用于Marangoni對流的研究.于藝紅等[11]通過激光紋影儀，直接在線監(jiān)測界面的湍動(dòng)對流結(jié)構(gòu)，記錄對流結(jié)構(gòu)的變化過程，為以后的研究提供堅(jiān)實(shí)的基礎(chǔ)；沙勇等[12,13]利用豎直流動(dòng)皂膜裝置，使用紋影光學(xué)方法觀察了由于丙酮從皂液中解吸，在微米級厚度皂膜上出現(xiàn)的滾筒狀Marangoni對流結(jié)構(gòu).

然而，如何解釋對流的微觀不穩(wěn)定性，以及宏觀的流動(dòng)演化過程，還需進(jìn)一步借助理論及數(shù)值計(jì)算的分析.例如，Straub[14]研究了失重條件下，表面張力在兩相流傳熱傳質(zhì)中的重要性，并理論分析了Marangoni效應(yīng)對傳熱的影響；沙勇等[15]結(jié)合流體動(dòng)力學(xué)條件，通過半經(jīng)驗(yàn)?zāi)Ｐ图袄碚摲治鲅芯苛薓arangoni對流的傳質(zhì)機(jī)理，并得出傳質(zhì)舍伍德數(shù)Sh與馬蘭戈尼數(shù)Ma的關(guān)聯(lián)；Hu等[16-18]采用潤滑分析法對干燥液滴進(jìn)行分析，并與流場的有限元算法相結(jié)合，得出表面活性劑能夠抑制Marangoni對流的結(jié)論.大量的實(shí)驗(yàn)與分析結(jié)果不僅證明了Marangoni現(xiàn)象改善了傳質(zhì)過程，而且還可以將實(shí)際傳遞速度與理論估計(jì)結(jié)果進(jìn)行比較，來量化臨界Ma數(shù)和傳遞速率的關(guān)系.在數(shù)值分析方面，早在1991年，唐澤眉等[19]采用有限元數(shù)值模擬，考慮邊界形狀的影響，得出流場、溫度場、表面壓強(qiáng)的分布；劉春元[20]建立熱毛細(xì)對流的數(shù)學(xué)模型并進(jìn)行數(shù)值模擬，分別得到有無相變時(shí)的速度、溫度分布；沙勇等[21]利用有機(jī)溶液能夠使水的表面張力降低的特性，人工誘發(fā)Marangoni對流，進(jìn)行數(shù)值模擬計(jì)算，研究了Marangoni對流對傳質(zhì)的影響；王貞濤等[22]建立了高溫基板上的液滴蒸發(fā)與高溫空氣中的懸浮液滴蒸發(fā)模型，研究蒸發(fā)過程中液滴內(nèi)部Marangoni對流非穩(wěn)態(tài)流動(dòng)模型；唐甜等[23]采用ALE方法建立固定接觸角蒸發(fā)模型，數(shù)值模擬研究了蒸發(fā)數(shù)Ec、接觸角θ對Marangoni對流的影響；Hu等[24]建立完善的蒸發(fā)模型，包括空氣中的蒸發(fā)擴(kuò)散，液體空氣界面中的蒸發(fā)冷卻，固體中的復(fù)合傳熱，浮力引起的對流和Marangoni對流，研究了加熱和不加熱基板上的蒸發(fā)特性.此外，還有很多學(xué)者對蒸發(fā)液滴的特性及誘導(dǎo)因素進(jìn)行了大量的仿真模擬[25-28]，分別獲得了在傾斜基板上、在油相中以及在激光照射下液滴內(nèi)Marangoni對流的流動(dòng)特性及其對傳熱傳質(zhì)過程的影響.

但是，目前的數(shù)值計(jì)算模型基板多為親水性材料，液滴接觸角在15 °～30 °之間，并沒有對疏水性液滴內(nèi)部的對流進(jìn)行深入研究.本文以固定基板上蒸發(fā)液滴為研究對象，基于COMSOL Multiphysics創(chuàng)建流體計(jì)算模型，通過與Abdullah等[29]的實(shí)驗(yàn)結(jié)果進(jìn)行對比驗(yàn)證分析，驗(yàn)證COMSOL Multiphysics數(shù)值模型的可靠性，并對接觸角θ從100 °到120 °變化的蒸發(fā)液滴進(jìn)行內(nèi)部穩(wěn)態(tài)流場特性分析，分析研究Marangoni對流的影響因素，為后續(xù)Marangoni對流的研究及優(yōu)化提供相應(yīng)的理論基礎(chǔ).

1 物理模型及網(wǎng)格劃分

1.1 物理模型

利用COMSOL軟件，建立的物理模型如圖1所示，主要分為液相區(qū)域、氣相區(qū)域兩部分.液相區(qū)域?yàn)樗?，水滴半徑R=0.7 mm，初始潤濕半徑為RC=0.6 mm，初始接觸角θ=2.09 rad；基底長度為1.5RC，氣相區(qū)域半徑為25R.液滴靜置于底層基板上，置于空氣中.氣、液兩相工質(zhì)分別為水、空氣.采用二維軸對稱模型，不考慮氣相物質(zhì)擴(kuò)散.

圖1 基底表面液滴蒸發(fā)示意圖

1.2 網(wǎng)格獨(dú)立性檢驗(yàn)

蒸發(fā)液滴計(jì)算域網(wǎng)格的劃分方法與模型的網(wǎng)格數(shù)量和網(wǎng)格質(zhì)量有著直接的聯(lián)系，而網(wǎng)格質(zhì)量的好壞對數(shù)值仿真預(yù)測的結(jié)果具有一定的影響[30].因此，本文通過考察接觸角θ=100 °，基板溫度為300.15K，t=2 s時(shí)液滴內(nèi)部最大流動(dòng)速度，對計(jì)算模型進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證.t=2 s時(shí)液滴內(nèi)部最大流動(dòng)速度結(jié)果見表1所示.

表1 網(wǎng)格無關(guān)性檢驗(yàn)研究

從表1可以看出，當(dāng)網(wǎng)格最小單元尺寸小于10-3mm時(shí)，對計(jì)算模型的結(jié)果影響較小，并且滿足穩(wěn)定性要求的步長.考慮計(jì)算時(shí)間和精度的影響，本文采取網(wǎng)格數(shù)量為5 962的網(wǎng)格.其中，網(wǎng)格劃分時(shí)，為了使界面處具有較好的網(wǎng)格質(zhì)量，主要對氣液界面和液固界面進(jìn)行網(wǎng)格加密處理.液滴內(nèi)部的網(wǎng)格對結(jié)果的影響不大，相關(guān)網(wǎng)格加密示意圖見圖2所示.

因此，本研究中氣液界面網(wǎng)格最小單元尺寸為0.000 217 mm，最大單元尺寸為0.041 1 mm，液固界面網(wǎng)格最小單元尺寸為0.000 749 mm，最大單元尺寸為0.026 2 mm，邊界層層數(shù)3層，邊界層拉伸因子1.2，所有域平均單元質(zhì)量為0.903，網(wǎng)格數(shù)量約為5 962.

(a)網(wǎng)格劃分示意圖 (b)界面網(wǎng)格加密示意圖 圖2 網(wǎng)格劃分及加密示意圖

2 數(shù)學(xué)模型及邊界條件

2.1 液滴內(nèi)部流動(dòng)及傳熱數(shù)學(xué)模型

本文選定COMSOL中層流流體流動(dòng)、傳熱、Marangoni效應(yīng)模型進(jìn)行蒸發(fā)液滴仿真數(shù)值模擬計(jì)算.

液滴內(nèi)部的流動(dòng)及傳熱模型可以由以下方程組描述.

(1)連續(xù)方程

(1)

式(1)中：ρ為流體密度，kg/m3；V為流體的速度，m/s.

(2)動(dòng)量方程(N-S方程)

(2)

式(2)中：P為流體微元體上的壓力，Pa；μ為流體的動(dòng)力粘度，Pa·s.

(3)能量方程

(3)

式(3)中：T為基板的溫度，K；cp為流體的比熱容，J/(kg·K)；k為流體的熱導(dǎo)率，W/(m·K).

Marangoni效應(yīng)在數(shù)值求解過程中是通過在N-S方程中添加額外的應(yīng)力作用來考慮的，即由于在與表面相切的s方向上存在的表面張力梯度dσ/ds會(huì)對自由表面流動(dòng)形成作用，因此就必然存在一個(gè)與s方向相反的剪切應(yīng)力來保證表面的受力平衡.該剪切應(yīng)力可能會(huì)改變邊界條件，從而改變流動(dòng)以及作用于液滴上的力.而產(chǎn)生該應(yīng)力作用的表面張力梯度可由溫度梯度、濃度梯度等諸多因素決定.研究人員發(fā)現(xiàn)，對于溫度梯度產(chǎn)生的表面張力變化，在大部分溫度范圍內(nèi)，表面張力隨溫度線性減小，直到臨界點(diǎn)時(shí)減小為零.因此，對于任何給定的流體，對其熱物理特性的dσ/dT進(jìn)行控制總是確定的，并且接近常數(shù)[31].

空氣中液滴表面的表面張力分布可以表示為：

(4)

假設(shè)熱交換主要由熱傳導(dǎo)控制，則液滴表面無熱量交換，那么利用拉普拉斯方程，可以得到：

(5)

式(4)、(5)中：R為液滴半徑，m；σ為液滴的表面張力系數(shù)，N/m；θ為接觸角，°；T為溫度，K；x1是x方向上任意一點(diǎn)的位置.

為使問題簡化，假設(shè)液滴是球形表面，法向方向上表面張力的差異非常小，法向的應(yīng)力完全由表面張力決定.

根據(jù)這些假設(shè)，球形液滴的切應(yīng)力邊界條件就變?yōu)椋?/p>

(6)

式(6)中：ρL為密度，kg/m3；vL為運(yùn)動(dòng)粘度，m2/s；u為速度，m/s.

那么根據(jù)方程得到作用于液滴上的力為：

(7)

除了法向阻力(第一項(xiàng))外，在沿液滴表面張力減小的方向上還有一個(gè)Marangoni力2πR2(dσ/dx1)作用在液滴上.因此，表面張力隨溫度增加而降低.均勻溫度梯度dT/dx1的存在就會(huì)使液滴在沿流體變熱的方向上承受一個(gè)大小為2πR2(-dσ/dT)(dT/dx1)的附加力.

2.2 邊界條件

(1)固液界面

固液兩相接觸界面采用無滑移邊界條件，溫度為常量.

V*=0

(8)

Θ=1

(9)

式(8)、(9)中：V*為壁面速度，Θ為熱通量.

(2)自由表面

切向力的平衡以及垂直于表面方向上熱通量的連續(xù)性，可用式(10)、(11)來表示：

(10)

(11)

式(10)、(11)中：Bi為畢渥數(shù)，Ma為馬蘭戈尼數(shù)，Pr為普朗特?cái)?shù).

依據(jù)上述建立的模型，對接觸角為100 °、110 °、120 °，及基板溫度為300.15 K、310.15 K、320.15 K時(shí)，蒸發(fā)液滴內(nèi)部的Marangoni對流流場分布及基板溫度對Marangoni效應(yīng)的影響進(jìn)行分析研究.

3 模擬結(jié)果及分析

3.1 模型驗(yàn)證

利用上述建立的模型，首先對Abdullah[29]等人的固著蒸發(fā)液滴的實(shí)驗(yàn)研究進(jìn)行同條件下的數(shù)值模擬，驗(yàn)證模型正確性.圖3為模擬得到的液滴內(nèi)部各點(diǎn)速度大小與原文實(shí)驗(yàn)結(jié)果的對比.由圖3可以看出，模型的計(jì)算結(jié)果與Abdullah[29]等的實(shí)驗(yàn)結(jié)果吻合較好，模型正確可靠.

圖3 液滴內(nèi)部速度的模擬與實(shí)驗(yàn)結(jié)果對比驗(yàn)證

3.2 蒸發(fā)液滴內(nèi)部Marangoni對流

3.2.1 Marangoni對流發(fā)展規(guī)律

初始狀態(tài)下，基板與液滴之間存在溫差，所以在基板和液滴之間會(huì)發(fā)生傳熱，在傳熱過程中產(chǎn)生溫度梯度，從而驅(qū)動(dòng)液滴界面的表面張力梯度，誘發(fā)Marangoni對流.圖4～8分別為t=0.05 s、0.1 s、0.15 s、0.2 s、0.3 s時(shí)液滴內(nèi)部的流場分布圖(a)及溫度分布圖(b).

從瞬態(tài)速度圖可以看出，初始階段(t=0～0.05 s時(shí))在液滴中下部出現(xiàn)小渦流，此時(shí)液滴邊緣界面處速度最大，且隨著時(shí)間變化界面速度不斷增大，渦流不斷變大.t=0～0.05 s內(nèi)，速度最大為4×10-6m/s，在t=0.2 s時(shí)達(dá)到1.2×10-5m/s，說明對流產(chǎn)生的機(jī)理是由于表面張力梯度的存在，在表面誘發(fā)，傳遞到液滴內(nèi)部.

從溫度分布圖可以看出，t=0.05 s時(shí)，熱量從基板傳遞到液滴內(nèi)部，三相界面處熱量傳遞速率比液滴內(nèi)部的傳遞速率大，這也是液滴中下部最先出現(xiàn)渦流的原因.t=0.2 s以后，液滴內(nèi)部速度最大區(qū)域始終在液滴中心及界面處；速度及溫度分布變化較小，說明t=0.2 s時(shí)內(nèi)部流動(dòng)達(dá)到穩(wěn)態(tài).

(a)速度矢量分布 (b)溫度分布圖4 t=0.05 s時(shí)液滴內(nèi)部流場分布及溫度分布

(a)速度矢量分布 (b)溫度分布圖5 t=0.1 s時(shí)液滴內(nèi)部流場分布及溫度分布

(a)速度矢量分布 (b)溫度分布圖6 t=0.15 s時(shí)液滴內(nèi)部流場分布及溫度分布

(a)速度矢量分布 (b)溫度分布圖7 t=0.2 s時(shí)液滴內(nèi)部流場分布及溫度分布

(a)速度矢量分布 (b)溫度分布圖8 t=0.3 s時(shí)液滴內(nèi)部流場分布及溫度分布

3.2.2 基板溫度對Marangoni對流的影響

圖9～11展示的是接觸角θ1=100 °和基板溫度分別為300.15 K、310.15 K、320.15 K時(shí)，液滴內(nèi)部Marangoni對流達(dá)到穩(wěn)態(tài)時(shí)的溫度分布及速度分布.從圖中可以看出，基板溫度越高，導(dǎo)致基板和液滴之間溫度梯度增大，所以對流速度越大，Marangoni效應(yīng)越明顯.在接觸角相同的情況下，由于基板溫度不同，液滴蒸發(fā)速度不同，基板溫度從300.15 K至320.15 K變化，最大速度由4.77×10-5m/s變?yōu)?.06×10-4m/s.

由圖可以很明顯地看出，軸對稱模型左側(cè)部分速度較大，即液滴中心速度比邊緣速度大，這是由于達(dá)到穩(wěn)態(tài)后，熱量向液滴中心傳遞速度比向液滴邊緣傳遞速度快.而相對應(yīng)的溫度分布圖中，隨著基板溫度增大，等溫線趨勢不變，但液滴內(nèi)最大溫度由299.63 K增至318.13 K，同水平線上液滴中心的溫度高于界面的溫度.

(a)溫度分布 (b)速度矢量分布圖9 T1=300.15 K時(shí)溫度分布及速度矢量分布

(a)溫度分布 (b)速度矢量分布圖10 T1=310.15 K時(shí)溫度分布及速度矢量分布

(a)溫度分布 (b)速度矢量分布圖11 T1=320.15 K時(shí)溫度分布及速度矢量分布

3.2.3 接觸角對Marangoni對流的影響

圖12～14分別展示的是接觸角為θ1=100 °、θ1=110 °、θ1=120 °和基板溫度為300.15 K時(shí)，液滴達(dá)到穩(wěn)態(tài)時(shí)的三維速度分布及速度矢量分布.從圖12～14可以看出，隨著接觸角的增大，傳熱區(qū)域也隨之增大，液滴內(nèi)部流動(dòng)速度相對增大，液滴內(nèi)部最大速度由4.77×10-5m/s增加至7×10-5m/s，可以發(fā)現(xiàn)接觸角對流動(dòng)速度有增強(qiáng)作用.

(a)三維速度分布 (b)速度矢量分布圖12 θ1=100 °三維速度分布及速度矢量分布

(a)三維速度分布 (b)速度矢量分布圖13 θ1=110 °三維速度分布及速度矢量分布

(a)三維速度分布 (b)速度矢量分布圖14 θ1=120 °三維速度分布及速度矢量分布

3.3 最大對流速度

本文分別對接觸角為100 °、110 °、120 °，基板溫度為300.15 K、310.15 K、320.15 K時(shí)的液滴進(jìn)行模擬，得出速度分布及溫度分布.初始階段液滴邊緣界面處速度最大，隨著時(shí)間變化，液滴達(dá)到穩(wěn)態(tài)，液滴內(nèi)部速度不斷增大，最終液滴內(nèi)部速度較大區(qū)域分布在液滴中心及界面處.

圖15顯示的是不同接觸角及基板溫度對液滴內(nèi)最大速度的影響.液滴最大流動(dòng)速度隨基板溫度的升高及接觸角的增大而增大.接觸角越大，流動(dòng)速度隨基板溫度增加的越快，這是由于液滴內(nèi)部與界面溫差增大和接觸角增大共同作用的結(jié)果.基板溫度升高時(shí)，Marangoni對流強(qiáng)化作用越明顯.

圖15 液滴接觸角與基板溫度對最大對流速度的影響

4 結(jié)論

(1)作瞬態(tài)研究時(shí)，初始階段液滴邊緣速度最大，這是由于界面溫度梯度引發(fā)的Marangoni對流在表面誘發(fā)，向內(nèi)部傳遞，t=0.2 s后內(nèi)部流動(dòng)達(dá)到穩(wěn)態(tài)，液滴中心速度逐漸增大，直至與界面速度相同.

(2)達(dá)到穩(wěn)態(tài)后，基板溫度固定不變時(shí)，接觸角從100 °到120 °變化，液滴內(nèi)部速度不斷變大，這是由于接觸角的增大使液相區(qū)域增大，傳熱區(qū)域增大，導(dǎo)致渦流變大，速度增強(qiáng)，說明在一定范圍內(nèi)接觸角的增大可強(qiáng)化Marangoni對流.

(3)達(dá)到穩(wěn)態(tài)后，接觸角固定不變時(shí)，基板溫度從300.15 K到320.15 K變化，與環(huán)境溫度形成的溫差變大，由溫度梯度引起的Marangoni效應(yīng)增強(qiáng)，液滴內(nèi)部流動(dòng)速度增大.基板溫度越高，Marangoni效應(yīng)形成的渦流結(jié)構(gòu)就越明顯，說明基板溫度的提高有利于相際間的傳熱.