时空调控微柱表面浸润性强化单气泡沸腾换热

陈宏霞，李林涵，王逸然，郭宇翔，刘霖

（华北电力大学能源动力与机械工程学院，北京102206）

引 言

核态沸腾是一种传热温差小、传热系数高的高效传热方式，在实际相变换热器中广泛存在。伴随电子芯片、燃料电池、航空航天以及核能等高能产业的发展，更多研究者开始关注如何调控、强化核态沸腾换热过程，从而进一步提高表面换热时的热通量。加工表面微结构及改变表面浸润性可直接影响核化气泡的动力学过程，是近年来热门的强化表面换热方法。

在换热器表面沉积纳米颗粒改变表面浸润性，是近年来应用广泛的方法。Deb等[1]和Hsu等[2]利用纳米SiO2颗粒涂层制备了从超亲水性到超疏水性变化的铜表面，发现涂层厚度会导致表面湿润性变化，当表面接触角θ接近150°时，气泡会合并形成气膜使传热恶化。Phan等[3]认为沸腾表面的最佳传热系数可以同时在θ接近90°或接近0°的情况下获得。Kumar等[4]将纳米颗粒厚度从180 nm增加到260 nm，则接触角从48°减小到42°，后者的换热性能更优。Quan等[5]观察发现适度亲水的纳米粒子被吸附在气泡界面上有效防止了气泡聚结，沉积的颗粒使表面粗糙度改善并使成核腔增多。除实验研究外，利用数值模拟研究池沸腾表面浸润性的方法也应用广泛。Zhao等[6]基于VOF方法建立二维亲疏水混合表面模型，研究了表面温度和传热系数以及气泡运动的变化规律。Zhang等[7]则设计了浸润性随温度改变的表面，低过热度时θ大有利于气泡成核，传热增强，高过热度时转为亲水性可提高CHF（临界热通量）。

纳米颗粒沉积可制备得到浸润性不同的表面，进一步与经过合理设计后的表面微结构耦合可使传热更加高效[8-12]。Može等[13]通过激光诱导与化学气相沉积制造了超疏水和超亲水微腔表面，发现具有更小直径微腔的超亲水表面传热性能更好。Zhang等[14]研究了润湿性和微结构的协同效应，设计了底部和侧壁表面θ均为113°，顶部θ分别为89°、9°和156°的微通道，发现疏水性降低气液相变能垒，增加成核位点；亲水性增强微结构的泵送能力，吸回更多液体到表面提高了CHF。Zhang等[15]则发现表面的孔隙活性会因蒸发量增大而得到改善。可见近年的研究已通过核态沸腾实验和数值模拟方法，清晰揭示了表面浸润性与微结构强化核态沸腾的机理。亲水性表面可促进气泡脱离，脱离直径减小，强化效果显著。疏水表面可促进气泡的成核及生长，低热通量下有可能具有强化效果；高热通量下多气泡容易合并成蒸汽膜恶化表面传热。因此可合理推断在核态沸腾的生长和脱离阶段若可实现浸润性的转换可获得更优的强化效果。陈彦君[16]基于Al2O3纳米颗粒在电场的电泳原理，研究高压电场（0～12 kV）控制流体中纳米颗粒的运动,认为纳米流体等效热导率随电场强度的增大而提高。Asadzadeh等[17]发现浓度为0.02%的纳米流体在12.5 kV电压下可获得最大自然对流传热系数。黄浩[18]则利用300 V/cm脉冲电场使纳米流体冷端温度发生周期性变化，同时研究了电场对换热表面纳米颗粒沉积的影响，认为可通过调控电场强度和时间获得合适的电泳沉积层，为实现核态沸腾过程中浸润性人为调控提供了一种新的可能。

同时，混杂浸润性表面的研究也日渐兴起[19-20]。Jo等[21]发现具有疏水点的亲水表面的BHT（沸腾传热系数）和CHF分别受疏水面积占总面积的比值r以及疏水点的直径、间距和数量的影响，并认为减小r的同时增大疏水点数量可保持高CHF和高BHT。Motezakker等[22]则制备了θ=20°和165°混杂的表面，认为r=38.5%的表面传热性能最佳。Betz等[23]设计了超疏水和超亲水区域并置的超双亲表面。Chen等[24]制作了具有疏水阵列的亲水表面，成核位点密度更高，气泡脱离直径减小为1/3，CHF提高90%。也有学者对表面微结构某一部分进行了浸润性修饰的研究，Jo等[25]利用SAM薄膜对微柱柱顶进行疏水性修饰，发现气泡更易在柱顶成核。Li等[26]则利用数值模拟方法设计了对柱顶进行浸润性修改的表面，发现柱底边缘和柱顶都易促使气泡成核，且气泡聚结与分离互不干扰，表面结构和混合润湿性的联合效应增强，核态沸腾传热得以强化。

近年来人们对微结构耦合均匀浸润性以及混杂浸润性等强化核态沸腾相变换热的研究较多，证明混杂浸润性表面可耦合亲、疏水表面各自对沸腾气泡的强化作用，进一步提高换热表面的换热性能。本文基于上述亲、疏水表面对核态沸腾的强化机理，提出分别从空间上、时间上对壁面浸润性进行调控的强化方法，利用CFD数值模拟研究三维微柱表面单气泡的核态沸腾过程；空间调控是保持柱顶以及整个表面始终更加亲水，时间调控是在某一时刻进行干预，使整个表面更加亲水，为提高换热表面对复杂换热环境的可控性和灵活性提供理论依据。

1 模型及方法

1.1 计算区域边界条件及网格划分

计算域几何尺寸与边界条件如图1(a)所示，为方便计算，仅对一半对称结构进行建模。硅板基底厚为50μm，单个微柱长宽高分别为50、50和150 μm，微柱间隙为100μm，成阵列分布；流体区域的高度为1500μm，设置平面y=0作为对称轴，计算域底面面积为1350μm×675μm，满足沸腾计算要求。基底下表面视为热通量58 kW/m2的加热表面，设置为第二类边界条件；y=0平面设置为轴对称边界条件。

计算域侧壁面边界条件的设置与计算沸腾单气泡的步骤有关：首先进行沸腾壁面纯导热计算，根据实验结果[27]计算成核点过热度达4 K时结束，获得成核壁面过热度及计算域内温度场和压力场；此阶段计算域侧壁面设置为绝热边界条件，上表面为压力出口，初始壁面温度为饱和温度。第二步将纯导热计算得到的流场作为沸腾过程初始流场，此时将三个侧壁面改为压力出口条件，并将侧壁面压力温度设置为对应压力出口条件下的压力和回流温度。最后在给定的成核点处人为放置种子气泡，在压力出口边界条件下开始沸腾过程计算至气泡完全脱离壁面。

网格划分情况如图1(b)所示，计算域均采用六面体网格划分，对于非气泡和远离壁面区域，采用较粗网格进行划分；水平方向网格长度为10μm，高度方向网格长度从10μm逐渐增加到17.28μm；对于整个基底区域和气泡整体区域进行网格加密，最小网格长度为2.5μm。

图1 计算域的设置Fig.1 Setting of computational domain

1.2 物性参数及工况设置

本文选择水作为沸腾工质，假设所有有关物性参数不随温度变化；沸腾过程涉及液态水和水蒸气两相，均假设为饱和温度和大气压力条件下的常物性不可压缩流体；整个沸腾换热过程为饱和沸腾，表面张力取为常值。涉及的材料与工质的物性参数如表1所示。

表1 数值模拟涉及工质的物性参数Table 1 Physical parameters employed in the numerical simulations

所涉及的沸腾模拟过程，具体涉及三种调控方式。空白试样为浸润性无调控的均匀浸润表面，壁面接触角始终为θ=48°。第一种调控是局部浸润，只在微柱顶部调控浸润性，保持柱顶壁面浸润性θ=20°，其他壁面浸润性θ=48°；第二种调控为梯度浸润，保持底部浸润性不变θ=48°，从侧壁到顶部壁面浸润性为梯度渐变，具体值为沿微柱侧壁高度每增加50μm减小至20°和10°，到柱顶为5°，该调控方式的表面亲水性最强；第一种和第二种浸润调控均为空间调控（图2）。第三种调控为时间调控，在t=0.40 ms时刻将表面初始浸润性θ=48°全部改变为θ=20°。选择水作为沸腾工质，假设所有有关物性参数不随温度变化；沸腾过程涉及液态水和水蒸气两相，均假设为饱和温度和大气压力条件下的常物性不可压缩流体；整个沸腾换热过程为饱和沸腾，表面张力取为常值。

图2 空间调控的两种方式Fig.2 Two cases of space modulation

1.3 数学模型及求解方法

本文基于CFD方法，使用商业软件Fluent15.0对不同浸润性调控下的单气泡核态沸腾过程进行三维数值模拟。采用已证明的数值计算模型[28]，除VOF方程和基本的控制方程外，气液界面蒸发模型、微层蒸发模型和表面张力模型以源项的形式加入动量方程。模型的验证方法是将光滑硅表面上的气泡形态演变与实验进行比较，以及将微柱结构表面上的气泡脱离直径与测量值进行比较，数值计算结果的准确性在可接受的误差范围内得到了验证[29]。此外，在与Utaka等[30]的实验相同的条件下，气泡单周期内微层蒸发量与蒸发总量的比例为39%，接近Utaka等值的36%，也验证了模型的有效性。

为保证计算过程界面的稳定性，将气液界面传热传质量赋值于界面临近网格，微层区域传热传质量赋值于微层区域顶层网格。在给定初始微层厚度与气泡半径关系的基础上，利用VOF方法[31]求解控制方程，通过界面自适应加密实现对气液界面的精确跟踪。

数值计算过程中采用压力和速度耦合以求解PISO算法，压力插值离散方案采用体积力加权格式；动量方程和能量方程的扩散相采用中心差分格式，对流项采用二阶迎风格式；体积分数方程采用界面重构插值方案；时间离散过程采用一阶隐式格式。采用变时间步长的方法，在种子气泡生长阶段取10-8s，脱离阶段增大到5×10-7s，可在节省计算时间的同时保证迭代结果准确。

2 结果与讨论

2.1 气泡动力学规律

图3展示了气泡体积随时间的变化规律,气泡先加速变大后增速减缓直到脱离。与均匀浸润相比，局部浸润与梯度浸润表面上单气泡体积和生长脱离周期都明显改变；单气泡周期由1.45 ms提前至1.42 ms和1.20 ms，分别缩短2.1%和17.2%；完全脱离时的气泡体积分别为0.382 mm3和0.434 mm3，增大4.3%和18.4%；表明空间调控可以缩短气泡周期、增大气泡体积和蒸发吸热量。时间调控与均匀浸润表面的气泡蒸发吸热量接近，体积都约为0.36 mm3；而前者脱离微柱时刻提前至1.14 ms,缩短21.4%，说明在初始成核点过热度和微柱高度密度相同的情况下，进行时间调控的效果主要体现在可明显缩短单气泡生长脱离周期。

图3 气泡体积V bubble随时间t的变化Fig.3 Variation of bubble volume V bubble with time t

图4(a)展示了气泡与换热表面基底的接触面积Ab-s随时间的变化规律。如Chen等[32]所述气泡的单个换热周期可分为生长阶段和脱离阶段；气泡在微柱表面的生长先由基底向四周和微柱上方扩展爬升，再由微柱上方向下不断包覆微柱，微柱间隙处的液体因受挤压而外流，在生长阶段Ab-s的变化是一个先增大后减小又再次回升的过程。本文定义气泡与基底接触面积减小后不再反弹的时刻，为开始进入脱离阶段的临界时刻，即生长阶段和脱离阶段的分界点。均匀浸润、局部浸润和梯度浸润表面上的气泡分别在0.495、0.52和0.505 ms进入脱离阶段；进行时间调控的时刻处于生长阶段，Ab-s先短暂减小，之后伴随气泡的继续生长而再次回升，进入脱离阶段的时刻由0.495 ms提前到0.45 ms，显著提前了气泡进入脱离阶段的时刻。气泡进入脱离阶段后，与基底的接触面积逐渐收缩减小，在均匀浸润表面于1.18 ms离开基底；梯度浸润曲线在0.8 ms之前基本与均匀浸润曲线重合，但之后迅速减小于1.04 ms变为0，显示了壁面亲水性增强对气泡脱离的促进作用。

图4 气泡与壁面接触面积随时间的变化Fig.4 Variation of contact area with time

图4(b)展示了气泡与微柱侧壁面接触面积Ab-p随时间的变化规律。生长阶段气泡包覆微柱并向周围扩张，接触面积不断增加；在脱离阶段气泡沿微柱向上爬升而逐渐减小。时间调控后表面更加亲水，气泡周围饱和液体沿界面内流，在微柱表面形成一层薄液膜，增强了气液界面蒸发作用，减小了气泡与侧壁面接触面积。梯度浸润表面的气泡因整体上更强的亲水性，薄液膜作用更加显著，Ab-p处于较低水平。

2.2 浸润性引起的界面变化

图5为t=0.50 ms时刻均匀浸润与时间调控、局部浸润和梯度浸润表面的气泡在y=0截面上的温度分布云图对比。在生长过程中气泡整个曲面受力平衡，当气泡曲面所受向上浮力大于与壁面的黏附力时气泡进入脱离阶段。为更好地分析微柱结构气泡曲面受力，将气泡所受表面张力以微柱柱顶高度为界分为方向竖直向下的Fsa与竖直向上的Fsb，如图5(a)所示。并简化认为Fsa代表壁面吸附力的物理作用，Fsb可代表浮力的物理作用。气泡合力可表示为：

根据Young-Laplace公式，气泡弯液面所受的表面张力Fs与曲率半径R的关系为：

式中,γ为表面张力系数，N/m；κ为界面曲率。

对比图5(a)、(b)，发现时间调控表面的气泡上半部分曲率变化很小，Fsa几乎不变；而如图5(a)、(b)中红色圆圈所示，壁面亲水性提高后气泡位于两个柱顶之间的气液界面下沉，曲率κ增大曲率半径减小，因此沿着气泡下半界面方向向上的Fsb增大；最终合力Fs增大促使气泡向上脱离。

浸润性的增强导致气泡周围饱和流体沿界面的回流，形成新薄液膜区域。对比观察图5(a)、(b)可发现回流引入的薄液膜取代了部分气泡与微柱壁面直接接触形成的干烧区，使相变换热量增加，蒸发速率提高，气泡换热增强、带走热量增大导致黑色圆圈内微柱间隙内出现了部分低温区域。对比图5(a)、(c)的柱顶部分，可发现局部浸润表面的柱顶附近也存在明显的薄液膜，而图5(d)梯度浸润中微柱上半部分更强的薄液膜换热使得蒸发速率有更大程度的提高，气泡相变换热量更大，导致圆圈内的低温区域更加明显。此外，薄液膜占据干烧区域，减小了气泡与微柱的直接接触面积，亲水性越强的表面上气泡与侧壁面接触面积越小。

图5 温度分布云图(t=0.50 ms，y=0截面)Fig.5 Temperature distribution with y=0 at t=0.5 ms

2.3 成核点过热度的演变规律

图5中展示了气泡在微柱结构上的受力以及浸润性变化带来的界面改变，也显示了气泡内部及周围流体域的温度分布情况，可知不同浸润性调控表面上气泡周围的热边界层厚度变化较小，可忽略；主要区别是气泡生长脱离过程中底部壁面过热度的变化。图6为成核点过热度随时间的变化曲线，对均匀浸润壁面，经历了生长阶段最初的波动之后，过热度随微层蒸发面积扩大、吸热增强而逐渐降低；在0.495 ms气泡进入脱离阶段后，底面接触面积逐渐减小，吸热减少导致热量在壁面成核点附近重新聚集，过热度开始回升直至下一个循环。局部浸润的过热度曲线相比均匀浸润的情况始终略低，过热度曲线最低的是梯度浸润表面，表明亲水性更强的壁面冷却效果更好；对于时间调控，亲水性增强使气泡加速脱离微柱，成核点处过热度回升较快，气泡体积和蒸发换热量与均匀浸润表面相差不大，所以脱离时刻的过热度也几乎相同。

图6 成核点过热度ΔT随时间t的变化Fig.6 Variation of superheatΔT with time t at the nucleation sites

2.4 气液界面蒸发速率和微层蒸发速率

图7(a)展示了气液界面蒸发速率qlv随时间的变化规律。在生长阶段随气泡的长大，气液界面蒸发速率先增大后减小，这是由于伴随气泡长大气液界面面积增大的同时越来越多的界面面积向热边界层的外缘移动，使得界面处的相变过热度及相变量减小；可见在气泡体积达到最大值前，气液界面处的换热速度已开始逐渐减小，并一直持续减小至气泡完全脱离微柱。对局部浸润，qlv在气泡爬升至亲水性柱顶之后便始终高于均匀浸润表面，梯度浸润更加明显，在0.205 ms时的最大速率比前两种情况提高约40%；时间调控使蒸发速率亦出现小幅度反弹，在0.03 ms内提高了10.5%，表明因亲水性增强导致周围流体内流而产生的薄液膜区域对气液界面蒸发速率的显著增强作用。

图7(b)展示了微层蒸发速率qml随时间的变化规律。曲线在达到最大值后，由于气泡沿微柱侧壁向上生长继而向下包覆后会再次在基底铺展，微层蒸发速率呈现先减小再回升的规律。然而，微层蒸发速率并不如图4(a)中底面接触面积在0.25 ms之后的回升幅度明显，原因是气泡在基底重新铺展的过程中，基底过热度如图6所示已经有所降低，微层蒸发的回升趋势受到限制。时间调控改变浸润性时气泡仍处于继续生长的阶段，qml也相应出现先减小后有所回升的波动；待气泡进入脱离阶段，微层换热持续减小。由于较强亲水性促进气泡脱离，时间调控和梯度浸润表面的微层蒸发速率曲线减小更快。

图7 气泡蒸发速率随时间t的变化Fig.7 Variation of evaporation rate with time t

在生长阶段微层蒸发起到重要作用，在脱离阶段气液界面速率逐渐占据主导；在0.70 ms时均匀浸润表面的气液界面速率为总速率的65.1%，局部浸润为68.4%，梯度浸润以及时间调控的占比则分别提升87%和85.5%。表面亲水性的增强使得气液界面速率占比增大，气泡在脱离阶段的换热能力有所增强。

2.5 单周期内换热性能的强化

综合考虑气泡体积即换热量以及时间的影响，定义从沸腾起始到某一时刻内的平均热流为α(t)

式中，Qml和Qlv分别为四种情况所对应的微层蒸发和气液界面蒸发换热量。

为了便于比较浸润性调控表面的气泡生长及脱离进程，将脱离时间进行归一化处理为t/td（td为四种表面各自的最终脱离时刻）。图8为热流α(t)随归一化时间t/td的变化规律；当t/td=1时，为气泡在一个单位周期内的当量热流。均匀浸润表面上的当量热流为0.262 W，在局部浸润、梯度浸润与时间浸润表面上则分别为0.278、0.374和0.323 W，分别增大6.1%、42.7%和23.3%，表明空间调控与时间调控都可显著强化单气泡周期内的沸腾换热。

图8 热流α(t)随t/t d的变化Fig.8 Variation of heat flowα(t)with t/t d

3 结 论

本文针对单气泡在150μm高度微柱表面的核态沸腾过程，采用CFD-VOF方法建立数值计算模型，对表面浸润性进行时间调控与空间调控，通过分析气泡的动力学与传热性能变化，得到如下结论。

（1）增强壁面的亲水性使得气泡下半部分界面曲率增大，气泡所受向上合力增强，促进了气泡脱离，显著缩短单气泡周期。相比均匀浸润表面的单气泡周期（1.45 ms），局部浸润表面、梯度浸润表面以及时间调控浸润表面单气泡周期为1.42、1.20和1.14 ms，分别缩短了2.1%和17.2%和21.4%。

（2）均匀浸润表面气泡生长阶段和脱离阶段内的平均热流分别为0.49 W和0.14 W，表明单气泡核态沸腾换热在生长阶段具有更高的换热效率；通过时空调控可增大生长阶段在单气泡换热周期内的时间占比，局部浸润、梯度浸润和时间调控的生长阶段用时占比分别延长了7.6%、23.8%和16.2%，单气泡周期内当量热流分别提升了6.1%、42.7%和23.3%，具有显著的强化效果。

（3）对比三种调控强化方法，梯度浸润表面具有最优的强化效果；但考虑实际微结构尺度下的梯度浸润实现难度较大，而基于目前电场对纳米颗粒的控制研究，时间调控浸润表面具有客观的强化效果和一定的可实现性。

符号说明

Ab-p——气泡与微柱侧壁及柱顶的接触面积，mm2

Ab-s——气泡与除微柱外的底面的接触面积，mm2

Fsa——表面吸附力，N

Fsb——浮力，N

q——蒸发速率，W

td——脱离阶段结束的时刻，ms

tg——开始脱离阶段的时刻，ms

αd——生长阶段内的热流，W

αg——脱离阶段内的热流，W

γ——表面张力系数，N/m

θ——表面的静态接触角,(°)

κ——气液界面的曲率，m-1

下角标

d——脱离阶段

g——生长阶段

lv——气液界面

ml——微层