基底厚度对蒸发液滴表面温度分布的影响

张凯，王依霖，徐学锋



基底厚度对蒸发液滴表面温度分布的影响

张凯，王依霖，徐学锋

（北京林业大学工学院，北京 100083）

蒸发液滴的表面温度分布对液滴的液体流动和颗粒沉积有着重要的影响。获得液滴表面温度目前主要采用数值计算方法。针对有限厚度基底上的蒸发液滴，分析了网格划分对液滴表面温度计算结果的影响。结果表明，相比于液滴边缘附近区域，液滴中心区域网格的细化对计算结果影响不大；而在接触线附近，相比于网格尺寸，网格细化区域大小对计算结果的影响也很小。利用数值方法研究了基底厚度对蒸发液滴表面温度分布特性的影响，发现随着基底厚度的改变液滴表面出现3种温度分布模式：（1）从液滴顶点到边缘处表面温度逐渐升高；（2）液滴表面温度非单调变化；（3）从液滴顶点到边缘处表面温度逐渐降低。考虑热传导路径长度和蒸发制冷的共同作用，对不同表面温度分布模式进行了解释，并获得了（R，）坐标平面上的表面温度分布模式相图。本文结果将有助于对液滴蒸发过程的理解，并为蒸发诱导自组装、喷墨印刷等技术提供理论依据。

液滴；蒸发；传热；表面温度分布；数值模拟

引 言

当今的工业生产和科学研究中，喷墨印刷[1-2]、喷射冷却[3]、自组装[4]以及薄膜凃层[5]等高新技术得到了广泛应用。在这些技术中，预测和控制蒸发液滴中颗粒的沉积过程和蒸发结束后颗粒的沉积形状是一个关键的问题。液滴蒸发过程中表面温度的不一致会引起表面张力梯度，进而改变液滴内部流场分布，并对液滴的颗粒沉积模式产生重要影响[6-14]。因此，蒸发液滴的表面温度分布已成为液滴蒸发中的一个重要研究内容，并得到国内外诸多专家学者的关注[7-17]。

蒸发过程中，液滴表面液体蒸发需消耗热量，这些热量由基底经由液滴内部传到液滴表面。Deegan等[6,18-19]认为，由于相对液滴表面的其他位置来说液滴顶点离基底的距离最长，其热传导路径也最长，因此液滴顶点的温度要低于液滴边缘处。Steinchen等[7]则认为液滴表面温度由蒸发强度决定。液滴表面蒸发强度由液滴顶部向液滴边缘不断增加，液滴边缘处由于蒸发强度最高因而温度最低。

Hu等[16,20-22]数值计算结果表明，蒸发液滴表面温度梯度方向与液滴接触角有关。当液滴接触角较大时，液滴表面温度由中心向边缘逐渐升高。当接触角减小到某一临界值以下时，液滴表面温度梯度方向会发生改变，表面温度由液滴中心向液滴边缘逐渐降低。Xu等[10-12 ]与Ristenpart等[13]进一步分析得出，液滴表面温度梯度方向发生改变时的临界接触角不仅取决于基底与液滴的相对热导率，而且与基底厚度和液滴接触线半径之比有关。此外，David等[8]和Dunn等[9]进行了基底热传递率对液滴温度分布影响方面的研究。

尽管针对蒸发液滴表面温度已进行了大量的研究，但是到目前为止，蒸发液滴表面温度分布依然是一个未被清楚了解的问题，而关于液滴表面温度变化规律众多学者依然有着不同的见解。其原因在于缺乏针对液滴表面温度变化规律的系统性研究。由于液滴蒸发问题的复杂性，目前普遍采用数值方法研究液滴表面温度问题。在数值计算中合理选择模型的网格划分方式，不仅可以减少计算工作量，还能尽可能地提高计算精度。本文系统研究了网格尺寸、网格细化区域大小以及网格细化位置对蒸发液滴表面温度计算结果的影响，其结论可用于指导进一步的研究和分析。进一步研究了基底相对厚度R对蒸发液滴表面温度分布的影响；同时针对不同的相对热导率R，得出了（R，）的相图。该研究有助于了解蒸发液滴中的Marangoni效应，并为控制液滴流场流动和液滴的颗粒沉积形式提供了依据。

1 蒸发液滴模型

如图1所示，在厚度为S、热导率为S的基底上，一个接触角为，接触线半径为、热导率为L的小液滴处于缓慢的蒸发中。由于液滴形状具有轴对称的特性，选择采用如图1所示的柱面坐标系（,）。

由于较小的毛细数和Bond数，液滴可以被近似地看作球冠形，因此液滴高度可表示为

对于缓慢蒸发的液滴，液滴表面的蒸发流量可近似为

式中，=1/2-/p，0与接触角、饱和蒸气压、相对湿度以及蒸汽扩散系数有关[6,19-20]。

在蒸发液滴中，对流传热相对于传导传热来说可以忽略，液滴内部热传导可以近似认为是一个准静态的过程[13,16,21,23]。因此，液滴内部温度满足拉普拉斯方程，用量纲1的形式表示为

式（3）的量纲1边界条件为：

为了分析网格划分以及基底厚度对于表面温度计算结果的影响，采用ANSYS软件对式（3）及边界条件式（4）～式（7）进行了数值求解。

2 网格划分对蒸发液滴表面温度分布的影响

通过研究网格细化位置对蒸发液滴表面温度分布的影响，发现液滴表面温度的计算精度和准确性主要与液滴边缘处的网格尺寸有关，而与液滴中心区域处网格尺寸关系不大。为了在保证计算精度的前提下减少计算工作量，应在液滴边缘处选取尽可能小的网格尺寸，而在液滴中心区域选择较大的网格尺寸。

接着研究了液滴边缘处对液滴表面温度分布的影响，发现液滴边缘处网格细化区域大小对计算结果的影响很小。因此，在液滴边缘处选取较小的细化区域，不仅能够保证计算的精度，而且可以减少计算的工作量。

最后研究网格尺寸对于蒸发液滴表面温度计算结果的影响。前面计算结果表明计算精度仅与液滴边缘处网格尺寸有关，因此，细化网格时仅选取液滴边缘处进行。选取了5种网格尺寸方案：① 全局网格尺寸为0.1；② 全局网格尺寸为0.01；③ 全局网格尺寸为0.01，且边缘处细化1次；④ 全局网格尺寸为0.01，且边缘处细化2次；⑤ 全局网格尺寸为0.01，且边缘细化3次。本节针对接触角=10º，相对热导率R=1，基底相对厚度R=0.2的蒸发液滴的表面温度进行了数值计算。

不同网格尺寸方案的液滴表面温度计算结果如图2所示。图中显示，液滴边缘处网格尺寸大小对计算结果的影响较大，但随着网格尺寸不断减小，相邻两种网格尺寸方案间计算结果的差异也不断变小。计算结果显示，液滴表面温度从液滴顶部到液滴边缘出现了先下降后上升的变化趋势，从而在液滴表面存在一个温度的最低值（图2）。随着网格尺寸的减小，液滴表面温度最低值不断下降，同时最低值出现的位置不断向接触线位置靠近。数值计算结果表明，网格单元尺寸为0.01，对接触线附近区域网格进行3次细化，即可满足如下收敛条件

3 基底厚度对蒸发液滴表面温度分布影响

基底相对厚度R对蒸发液滴表面温度分布有显著影响。Xu等[10]通过渐进分析发现，液滴温度梯度变化的临界角与基底相对厚度有关。为了进一步研究基底相对厚度对液滴表面温度分布的影响，本文针对接触角=10º，相对热导率R=0.5，基底相对厚度R分别为0、0.01、0.02、0.1、0.2、0.4、0.6、0.8、1的蒸发液滴的表面温度进行了数值计算。

3.1 蒸发液滴表面温度分布特性

如图3所示，当相对热传导率R、液滴接触角一定时，随着基底相对厚度R的变化，液滴表面温度分布趋势也会发生变化。当基底相对厚度R趋近于0时，液滴表面温度随着与液滴顶点距离的增大而逐渐升高。当基底相对厚度R逐渐增大时，液滴表面温度依次出现以下分布趋势：随着与液滴顶点距离的增大非单调变化；随着与液滴顶点距离的增大而逐渐下降；随着与液滴顶点距离的增大非单调变化；随着与液滴顶点距离的增大而逐渐升高。

当基底厚度较小（R≤0.6）时，液滴表面温度随基底厚度的变化趋势可以用蒸发制冷作用来解释。蒸发过程中，液体在液滴表面的蒸发需要消耗热量，从而引起液滴表面温度降低，而维持蒸发所需热量由基底经由液滴内部传导到液滴表面。假定在液滴和基底中热传导方向沿轴向上（即一维热传导假设[10]），因而液滴表面温度由蒸发制冷强度和热传导路径长度决定。

当R趋近于0时，液滴顶点与液滴边缘处的热传导路径差别较大，此时热传导路径长度对于液滴表面温度分布起决定作用。由于表面温度会随着热传导路径长度（即液滴表面高度）的增加而降低，故出现液滴表面温度随与液滴顶点距离增加而单调升高的分布。当R逐渐增大时，液滴顶点与边缘处热传导路径长度的差别逐渐变小，从而蒸发制冷对液滴表面温度分布的影响逐渐增加。由于液滴表面蒸发强度随与液滴顶点距离增加而增大，蒸发制冷作用决定的表面温度分布与热传导长度的影响相反，温度会随与液滴顶点距离的增加而降低。因此，随着R的增大，蒸发制冷和热传导路径的共同作用导致液滴表面温度发生非单调变化，并进而随着与液滴顶点距离的增大而降低（图3）。

但是，如图3所示，当基底厚度继续增大时，液滴表面温度分布又重新出现非单调变化，并逐渐发展为随着与液滴顶点距离的增大而升高。这一现象无法用上面的模型进行解释。出现这种现象的原因，是由于基底厚度较大时基底和液滴内部的温度场分布较为复杂，简单的一维热传导假设已经不再适用。对于较厚基底上蒸发液滴表面温度分布的研究，需要对本文的模型做进一步的改进。

3.2 蒸发液滴表面温度分布的相图（R，）

针对不同相对热传导率R，分别计算了液滴表面温度分布模式随基底相对厚度R和液滴接触角的变化，从而建立了坐标平面（R，）上表面温度分布模式的相图，如图4～图6所示。

由计算结果可知，当R一定时，随着基底相对厚度R和液滴接触角的改变，液滴表面温度分布模式会经历3种状态：① 在区域Ⅰ中，从液滴顶点到边缘区域，液滴表面温度会逐渐升高；② 在区域Ⅱ中，从液滴顶点到边缘区域，液滴表面温度会出现非单调的变化；③ 在区域Ⅲ中，从液滴顶点到边缘区域，液滴表面温度会逐渐下降。

此外，由图4～图6可以看出，随着R的增大，相图中各区域转化的临界角C都出现先上升后下降的变化。上升部分可用热传导路径长度的变化解释。接触角增大导致液滴表面各处热传导路径长度的差异增大，而基底厚度增大会导致传导路径长度的差异变小。当基底厚度较大时，会使得在较大的接触角时热传导路径长度的影响已经不再起决定作用，从而温度分布模式由区域1转变到区域2。临界角C随R的增大而下降部分，无法用目前的模型解释。这再次表明，当基底厚度较大时，采用的一维热传导模型就不再适用。

4 结 论

本文采用有限厚度基底的液滴蒸发模型，利用ANSYS对蒸发液滴表面温度分布进行了数值模拟，研究了网格细化位置、网格细化区域大小、网格尺寸对液滴表面温度计算的影响，并进一步研究了相对基底厚度R和接触角对液滴表面温度分布的影响，得出了如下结论。

（1）液滴边缘处网格尺寸对数值计算结果影响较大。为保证计算精度，建立模型时应对液滴边缘处网格进行细化，同时为节省计算量可选取较小的细化区域。

（2）随着相对基底厚度R的变化，液滴表面温度分布会出现3种模式：① 从液滴顶点到边缘表面温度逐渐升高；② 从液滴顶点到边缘表面温度非单调变化；③ 从液滴顶点到边缘表面温度逐渐下降。这3种状态可用热传导路径长度和蒸发制冷作用来解释。

（3）建立了基底相对厚度R和接触角坐标平面（R，）上液滴表面温度分布模式的相图，发现随着R的增大，相图中各区域转化的临界角C先上升后下降。

符 号 说 明

Bo——Bond数 Ca——毛细数 er——径向单位矢量 ez——轴向单位矢量 H——蒸发潜热 hR——基底相对厚度，hR=hS/R hS——基底厚度，mm J0——蒸发常数，与接触角、饱和蒸气压、相对湿度以及蒸汽扩散系数有关 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)；教育部新世纪优秀人才支持计划项目(NCET-12-0786)；高等学校博士学科点专项科研基金项目(20120014120017)。

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

联系人：徐学锋。第一作者：张凯（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).