張凱，王依霖，徐學鋒

?

基底厚度對蒸發(fā)液滴表面溫度分布的影響

張凱，王依霖，徐學鋒

（北京林業(yè)大學工學院，北京 100083）

蒸發(fā)液滴的表面溫度分布對液滴的液體流動和顆粒沉積有著重要的影響。獲得液滴表面溫度目前主要采用數(shù)值計算方法。針對有限厚度基底上的蒸發(fā)液滴，分析了網(wǎng)格劃分對液滴表面溫度計算結果的影響。結果表明，相比于液滴邊緣附近區(qū)域，液滴中心區(qū)域網(wǎng)格的細化對計算結果影響不大；而在接觸線附近，相比于網(wǎng)格尺寸，網(wǎng)格細化區(qū)域大小對計算結果的影響也很小。利用數(shù)值方法研究了基底厚度對蒸發(fā)液滴表面溫度分布特性的影響，發(fā)現(xiàn)隨著基底厚度的改變液滴表面出現(xiàn)3種溫度分布模式：（1）從液滴頂點到邊緣處表面溫度逐漸升高；（2）液滴表面溫度非單調(diào)變化；（3）從液滴頂點到邊緣處表面溫度逐漸降低。考慮熱傳導路徑長度和蒸發(fā)制冷的共同作用，對不同表面溫度分布模式進行了解釋，并獲得了（R，）坐標平面上的表面溫度分布模式相圖。本文結果將有助于對液滴蒸發(fā)過程的理解，并為蒸發(fā)誘導自組裝、噴墨印刷等技術提供理論依據(jù)。

液滴；蒸發(fā)；傳熱；表面溫度分布；數(shù)值模擬

引 言

當今的工業(yè)生產(chǎn)和科學研究中，噴墨印刷[1-2]、噴射冷卻[3]、自組裝[4]以及薄膜凃層[5]等高新技術得到了廣泛應用。在這些技術中，預測和控制蒸發(fā)液滴中顆粒的沉積過程和蒸發(fā)結束后顆粒的沉積形狀是一個關鍵的問題。液滴蒸發(fā)過程中表面溫度的不一致會引起表面張力梯度，進而改變液滴內(nèi)部流場分布，并對液滴的顆粒沉積模式產(chǎn)生重要影響[6-14]。因此，蒸發(fā)液滴的表面溫度分布已成為液滴蒸發(fā)中的一個重要研究內(nèi)容，并得到國內(nèi)外諸多專家學者的關注[7-17]。

蒸發(fā)過程中，液滴表面液體蒸發(fā)需消耗熱量，這些熱量由基底經(jīng)由液滴內(nèi)部傳到液滴表面。Deegan等[6,18-19]認為，由于相對液滴表面的其他位置來說液滴頂點離基底的距離最長，其熱傳導路徑也最長，因此液滴頂點的溫度要低于液滴邊緣處。Steinchen等[7]則認為液滴表面溫度由蒸發(fā)強度決定。液滴表面蒸發(fā)強度由液滴頂部向液滴邊緣不斷增加，液滴邊緣處由于蒸發(fā)強度最高因而溫度最低。

Hu等[16,20-22]數(shù)值計算結果表明，蒸發(fā)液滴表面溫度梯度方向與液滴接觸角有關。當液滴接觸角較大時，液滴表面溫度由中心向邊緣逐漸升高。當接觸角減小到某一臨界值以下時，液滴表面溫度梯度方向會發(fā)生改變，表面溫度由液滴中心向液滴邊緣逐漸降低。Xu等[10-12 ]與Ristenpart等[13]進一步分析得出，液滴表面溫度梯度方向發(fā)生改變時的臨界接觸角不僅取決于基底與液滴的相對熱導率，而且與基底厚度和液滴接觸線半徑之比有關。此外，David等[8]和Dunn等[9]進行了基底熱傳遞率對液滴溫度分布影響方面的研究。

盡管針對蒸發(fā)液滴表面溫度已進行了大量的研究，但是到目前為止，蒸發(fā)液滴表面溫度分布依然是一個未被清楚了解的問題，而關于液滴表面溫度變化規(guī)律眾多學者依然有著不同的見解。其原因在于缺乏針對液滴表面溫度變化規(guī)律的系統(tǒng)性研究。由于液滴蒸發(fā)問題的復雜性，目前普遍采用數(shù)值方法研究液滴表面溫度問題。在數(shù)值計算中合理選擇模型的網(wǎng)格劃分方式，不僅可以減少計算工作量，還能盡可能地提高計算精度。本文系統(tǒng)研究了網(wǎng)格尺寸、網(wǎng)格細化區(qū)域大小以及網(wǎng)格細化位置對蒸發(fā)液滴表面溫度計算結果的影響，其結論可用于指導進一步的研究和分析。進一步研究了基底相對厚度R對蒸發(fā)液滴表面溫度分布的影響；同時針對不同的相對熱導率R，得出了（R，）的相圖。該研究有助于了解蒸發(fā)液滴中的Marangoni效應，并為控制液滴流場流動和液滴的顆粒沉積形式提供了依據(jù)。

1 蒸發(fā)液滴模型

如圖1所示，在厚度為S、熱導率為S的基底上，一個接觸角為，接觸線半徑為、熱導率為L的小液滴處于緩慢的蒸發(fā)中。由于液滴形狀具有軸對稱的特性，選擇采用如圖1所示的柱面坐標系（,）。

由于較小的毛細數(shù)和Bond數(shù)，液滴可以被近似地看作球冠形，因此液滴高度可表示為

對于緩慢蒸發(fā)的液滴，液滴表面的蒸發(fā)流量可近似為

式中，=1/2-/p，0與接觸角、飽和蒸氣壓、相對濕度以及蒸汽擴散系數(shù)有關[6,19-20]。

在蒸發(fā)液滴中，對流傳熱相對于傳導傳熱來說可以忽略，液滴內(nèi)部熱傳導可以近似認為是一個準靜態(tài)的過程[13,16,21,23]。因此，液滴內(nèi)部溫度滿足拉普拉斯方程，用量綱1的形式表示為

式（3）的量綱1邊界條件為：

為了分析網(wǎng)格劃分以及基底厚度對于表面溫度計算結果的影響，采用ANSYS軟件對式（3）及邊界條件式（4）～式（7）進行了數(shù)值求解。

2 網(wǎng)格劃分對蒸發(fā)液滴表面溫度分布的影響

通過研究網(wǎng)格細化位置對蒸發(fā)液滴表面溫度分布的影響，發(fā)現(xiàn)液滴表面溫度的計算精度和準確性主要與液滴邊緣處的網(wǎng)格尺寸有關，而與液滴中心區(qū)域處網(wǎng)格尺寸關系不大。為了在保證計算精度的前提下減少計算工作量，應在液滴邊緣處選取盡可能小的網(wǎng)格尺寸，而在液滴中心區(qū)域選擇較大的網(wǎng)格尺寸。

接著研究了液滴邊緣處對液滴表面溫度分布的影響，發(fā)現(xiàn)液滴邊緣處網(wǎng)格細化區(qū)域大小對計算結果的影響很小。因此，在液滴邊緣處選取較小的細化區(qū)域，不僅能夠保證計算的精度，而且可以減少計算的工作量。

最后研究網(wǎng)格尺寸對于蒸發(fā)液滴表面溫度計算結果的影響。前面計算結果表明計算精度僅與液滴邊緣處網(wǎng)格尺寸有關，因此，細化網(wǎng)格時僅選取液滴邊緣處進行。選取了5種網(wǎng)格尺寸方案：① 全局網(wǎng)格尺寸為0.1；② 全局網(wǎng)格尺寸為0.01；③ 全局網(wǎng)格尺寸為0.01，且邊緣處細化1次；④ 全局網(wǎng)格尺寸為0.01，且邊緣處細化2次；⑤ 全局網(wǎng)格尺寸為0.01，且邊緣細化3次。本節(jié)針對接觸角=10o，相對熱導率R=1，基底相對厚度R=0.2的蒸發(fā)液滴的表面溫度進行了數(shù)值計算。

不同網(wǎng)格尺寸方案的液滴表面溫度計算結果如圖2所示。圖中顯示，液滴邊緣處網(wǎng)格尺寸大小對計算結果的影響較大，但隨著網(wǎng)格尺寸不斷減小，相鄰兩種網(wǎng)格尺寸方案間計算結果的差異也不斷變小。計算結果顯示，液滴表面溫度從液滴頂部到液滴邊緣出現(xiàn)了先下降后上升的變化趨勢，從而在液滴表面存在一個溫度的最低值（圖2）。隨著網(wǎng)格尺寸的減小，液滴表面溫度最低值不斷下降，同時最低值出現(xiàn)的位置不斷向接觸線位置靠近。數(shù)值計算結果表明，網(wǎng)格單元尺寸為0.01，對接觸線附近區(qū)域網(wǎng)格進行3次細化，即可滿足如下收斂條件

3 基底厚度對蒸發(fā)液滴表面溫度分布影響

基底相對厚度R對蒸發(fā)液滴表面溫度分布有顯著影響。Xu等[10]通過漸進分析發(fā)現(xiàn)，液滴溫度梯度變化的臨界角與基底相對厚度有關。為了進一步研究基底相對厚度對液滴表面溫度分布的影響，本文針對接觸角=10o，相對熱導率R=0.5，基底相對厚度R分別為0、0.01、0.02、0.1、0.2、0.4、0.6、0.8、1的蒸發(fā)液滴的表面溫度進行了數(shù)值計算。

3.1 蒸發(fā)液滴表面溫度分布特性

如圖3所示，當相對熱傳導率R、液滴接觸角一定時，隨著基底相對厚度R的變化，液滴表面溫度分布趨勢也會發(fā)生變化。當基底相對厚度R趨近于0時，液滴表面溫度隨著與液滴頂點距離的增大而逐漸升高。當基底相對厚度R逐漸增大時，液滴表面溫度依次出現(xiàn)以下分布趨勢：隨著與液滴頂點距離的增大非單調(diào)變化；隨著與液滴頂點距離的增大而逐漸下降；隨著與液滴頂點距離的增大非單調(diào)變化；隨著與液滴頂點距離的增大而逐漸升高。

當基底厚度較小（R≤0.6）時，液滴表面溫度隨基底厚度的變化趨勢可以用蒸發(fā)制冷作用來解釋。蒸發(fā)過程中，液體在液滴表面的蒸發(fā)需要消耗熱量，從而引起液滴表面溫度降低，而維持蒸發(fā)所需熱量由基底經(jīng)由液滴內(nèi)部傳導到液滴表面。假定在液滴和基底中熱傳導方向沿軸向上（即一維熱傳導假設[10]），因而液滴表面溫度由蒸發(fā)制冷強度和熱傳導路徑長度決定。

當R趨近于0時，液滴頂點與液滴邊緣處的熱傳導路徑差別較大，此時熱傳導路徑長度對于液滴表面溫度分布起決定作用。由于表面溫度會隨著熱傳導路徑長度（即液滴表面高度）的增加而降低，故出現(xiàn)液滴表面溫度隨與液滴頂點距離增加而單調(diào)升高的分布。當R逐漸增大時，液滴頂點與邊緣處熱傳導路徑長度的差別逐漸變小，從而蒸發(fā)制冷對液滴表面溫度分布的影響逐漸增加。由于液滴表面蒸發(fā)強度隨與液滴頂點距離增加而增大，蒸發(fā)制冷作用決定的表面溫度分布與熱傳導長度的影響相反，溫度會隨與液滴頂點距離的增加而降低。因此，隨著R的增大，蒸發(fā)制冷和熱傳導路徑的共同作用導致液滴表面溫度發(fā)生非單調(diào)變化，并進而隨著與液滴頂點距離的增大而降低（圖3）。

但是，如圖3所示，當基底厚度繼續(xù)增大時，液滴表面溫度分布又重新出現(xiàn)非單調(diào)變化，并逐漸發(fā)展為隨著與液滴頂點距離的增大而升高。這一現(xiàn)象無法用上面的模型進行解釋。出現(xiàn)這種現(xiàn)象的原因，是由于基底厚度較大時基底和液滴內(nèi)部的溫度場分布較為復雜，簡單的一維熱傳導假設已經(jīng)不再適用。對于較厚基底上蒸發(fā)液滴表面溫度分布的研究，需要對本文的模型做進一步的改進。

3.2 蒸發(fā)液滴表面溫度分布的相圖（R，）

針對不同相對熱傳導率R，分別計算了液滴表面溫度分布模式隨基底相對厚度R和液滴接觸角的變化，從而建立了坐標平面（R，）上表面溫度分布模式的相圖，如圖4～圖6所示。

由計算結果可知，當R一定時，隨著基底相對厚度R和液滴接觸角的改變，液滴表面溫度分布模式會經(jīng)歷3種狀態(tài)：① 在區(qū)域Ⅰ中，從液滴頂點到邊緣區(qū)域，液滴表面溫度會逐漸升高；② 在區(qū)域Ⅱ中，從液滴頂點到邊緣區(qū)域，液滴表面溫度會出現(xiàn)非單調(diào)的變化；③ 在區(qū)域Ⅲ中，從液滴頂點到邊緣區(qū)域，液滴表面溫度會逐漸下降。

此外，由圖4～圖6可以看出，隨著R的增大，相圖中各區(qū)域轉化的臨界角C都出現(xiàn)先上升后下降的變化。上升部分可用熱傳導路徑長度的變化解釋。接觸角增大導致液滴表面各處熱傳導路徑長度的差異增大，而基底厚度增大會導致傳導路徑長度的差異變小。當基底厚度較大時，會使得在較大的接觸角時熱傳導路徑長度的影響已經(jīng)不再起決定作用，從而溫度分布模式由區(qū)域1轉變到區(qū)域2。臨界角C隨R的增大而下降部分，無法用目前的模型解釋。這再次表明，當基底厚度較大時，采用的一維熱傳導模型就不再適用。

4 結 論

本文采用有限厚度基底的液滴蒸發(fā)模型，利用ANSYS對蒸發(fā)液滴表面溫度分布進行了數(shù)值模擬，研究了網(wǎng)格細化位置、網(wǎng)格細化區(qū)域大小、網(wǎng)格尺寸對液滴表面溫度計算的影響，并進一步研究了相對基底厚度R和接觸角對液滴表面溫度分布的影響，得出了如下結論。

（1）液滴邊緣處網(wǎng)格尺寸對數(shù)值計算結果影響較大。為保證計算精度，建立模型時應對液滴邊緣處網(wǎng)格進行細化，同時為節(jié)省計算量可選取較小的細化區(qū)域。

（2）隨著相對基底厚度R的變化，液滴表面溫度分布會出現(xiàn)3種模式：① 從液滴頂點到邊緣表面溫度逐漸升高；② 從液滴頂點到邊緣表面溫度非單調(diào)變化；③ 從液滴頂點到邊緣表面溫度逐漸下降。這3種狀態(tài)可用熱傳導路徑長度和蒸發(fā)制冷作用來解釋。

（3）建立了基底相對厚度R和接觸角坐標平面（R，）上液滴表面溫度分布模式的相圖，發(fā)現(xiàn)隨著R的增大，相圖中各區(qū)域轉化的臨界角C先上升后下降。

符 號 說 明

Bo——Bond數(shù) Ca——毛細數(shù) er——徑向單位矢量 ez——軸向單位矢量 H——蒸發(fā)潛熱 hR——基底相對厚度，hR=hS/R hS——基底厚度，mm J0——蒸發(fā)常數(shù)，與接觸角、飽和蒸氣壓、相對濕度以及蒸汽擴散系數(shù)有關 kL——液滴熱導率，cal·cm-1·s-1·K-1 kR——基底相對熱導率，kR=kS/kL kS——基底熱導率，cal·cm-1·s-1·K-1 R——接觸線半徑，mm ——量綱1化的液滴溫度 ——量綱1化的基底溫度 q——液滴接觸角

References

[1] Calvert P. Inkjet printing for materials and devices [J]., 2001, 13(10): 3299-3305

[2] Kim D, Jeong S, Park B K,. Direct writing of silver conductive patterns: improvement of film morphology and conductance by controlling solvent compositions [J]., 2006, 89: 264101

[3] Jia W, Qiu H H. Experimental investigation of droplet dynamics and heat transfer in spray cooling [J]., 2003, 27: 829-838

[4] Schnall-Levin M, Lauga E, Brenner M P. Self-assembly of spherical particles on an evaporating sessile droplet [J]., 2006, 22: 4547-4551

[5] Kimura M, Misner M J, Xu T,. Long-range ordering of diblock copolymers induced by droplet pinning [J]., 2003, 19(23): 9910-9913

[6] Deegan R D, Bakajin O, Dupont T F,. Contact line deposits in an evaporating drop [J]., 2000, 62: 756-765

[7] Steinchen A, Sefiane K J. Self-organised Marangoni motion at evaporating drops or in capillary menisci - thermohydrodynamical model [J]., 2005, 30(1): 39-51

[8] David S，Sefiane K, Tadrist L. Experimental investigation of the effect of thermal properties of the substrate in the wetting and evaporation of sessile drops [J].:, 2007, 298(1/2): 108-114

[9] Dunn G J, Wilson S K, Duffy B R,. A mathematical model for the evaporation of a thin sessile liquid droplet: comparison between experiment and theory [J].:, 2008, 323(1/2/3): 50-55

[10] Xu X F, Luo J B, Guo D. Criterion for reversal of thermal Marangoni flow in drying drops [J]., 2010, 26(3): 1918-1922

[11] Xu X F, Luo J B, Guo D. Radial-velocity profile along the surface of evaporating liquid droplets[J]., 2012, 8: 5797-5803

[12] Xu X F, Luo J B, Marangoni flow in an evaporating water droplet [J].., 2007, 91: 124102

[13] Ristenpart W D, Kim P G, Dmingues C,. Influence of substrate conductivity on circulation reversal in evaporating drops[J]., 2007, 99(23): 234502

[14] Lin Handan(林晗丹), Yu Guangxiong(余光雄), Jiao Fangjian(焦放健), Sha Yong(沙勇). Marangoni convection at droplet interface during gas-liquid mass transfer [J].(化工學報), 2013, 64(8): 2846-2852

[15] Hu Dinghua(胡定華), Wu Huiying(吳慧英), Wu Xinyu(吳信宇). Surface temperature distribution of droplet heated by local microheater [J].(化工學報), 2011, 62(11): 3039-3045

[16] Hu H, Larson R G. Analysis of the effects of Marangoni stresses on the microflow in an evaporating sessile droplet [J]., 2005, 21(9): 3972-3980

[17] Xin Juanjuan(辛娟娟), Zhou Zhifu(周致富), Xin Hui(辛慧), Chen Bin(陳斌), Wang Guoxiang(王國祥), Guo Liejin(郭烈錦). Validation analysis of different mass transfer formula in single droplet evaporation model [J].(化工學報), 2012, 63(6): 1704-1708

[18] Deegan R D. Pattern formation in drying drops [J]., 2000, 61: 475-485

[19] Deegan R D, Bakajin O, Dupont T F,. Capillary flow as the cause of ring stains from dried liquid drops [J]., 1997, 389: 827-829

[20] Hu H, Larson R G. Evaporation of a sessile droplet on a substrate [J]., 2002, 106(6): 1334-1344

[21] Hu H, Larson R G. Analysis of the microfluid flow in an evaporating sessile droplet [J]., 2005, 21(9): 3963-3971

[22] Hu H, Larson R G. Marangoni effect reverses coffee-ring depositions [J]., 2006, 110: 7090-7094

[23] Girard F, Antoni M, Sefiane K. On the effect of Marangoni flow on evaporation rates of heated water drops [J]., 2008, 24(17): 9207-9210

Influence of substrate thickness on temperature distribution along surface of drying droplets

ZHANG Kai，WANG Yilin，XU Xuefeng

(School of Technology，Beijing Forestry University，Beijing 100083，China)

The temperature distribution along the surface of drying droplet has significant influence on fluid flow and particle deposition of the droplet. Currently, the surface temperature distribution of the droplet is often studied by numerical methods. For drying droplets on finite thickness substrates, the effect of meshing on calculated surface temperature is discussed. Compared to the region near the droplet edge, mesh refinement in the central region of the droplet has negligible effect on numerical results. The influence of the size of mesh refinement region near the droplet edge on droplet surface temperature can also be neglected compared to the level of mesh refinement. By studying the effects of relative substrate thicknessRon droplet surface temperature, three different states with changing substrate thickness are found. (1) From center to edge, surface temperature increases monotonically. (2) From center to edge, surface temperature changes non-monotonically. (3) From center to edge, surface temperature decreases monotonically. The different patterns of droplet surface temperature can be explained by considering the combined effects of heat conduction path length and evaporative cooling. The surface temperature “phase diagram” on parameters(R,) with different relative thermal conductivities is obtained. The results will be helpful to understanding droplet evaporation and providing theoretical basis for evaporation induced self-assembly and ink jet printing.

droplet; evaporation; heat transfer; surface temperature distribution; numerical simulation

2014-05-20.

XU Xuefeng,associate professor, xuxuefeng @bjfu. edu.cn

10. 11949/j.issn.0438-1157.20140745

TK 124

A

0438—1157（2015）02—0703—06

國家自然科學基金項目(51275050)；教育部新世紀優(yōu)秀人才支持計劃項目(NCET-12-0786)；高等學校博士學科點專項科研基金項目(20120014120017)。

2014-05-20收到初稿，2014-09-18收到修改稿。

聯(lián)系人：徐學鋒。第一作者：張凱（1995—），男，本科。

supported by the National Natural Science Foundation of China (51275050), the Program for New Century Excellent Talents in University (NCET-12-0786) and the Specialized Research Fund for the Doctoral Program of Higher Education (20120014120017).