方肋微通道内流动沸腾的气泡动态与传热特性分析

申宇，潘振海，吴慧英

（上海交通大学机械与动力工程学院，上海200240）

随着近年来微电子技术和MEMS工艺在诸多工业领域的飞速发展，电子器件在尺寸缩小的同时，功耗却不断攀升。这导致电子设备表面的热通量激增，严重影响系统性能甚至直接烧毁电子器件[1]。带肋微通道热沉的流动沸腾换热在高热通量芯片热控制方面具有显著优势，得到了广泛关注[2-10]：与常规通道不同，带肋微通道的换热比表面积大，结构紧凑，便于封装；流动沸腾利用液体相变潜热，换热效率高，并能维持较好的壁面均温性。此外，相较于平滑微通道，其内部的微肋结构还具有强化流体扰动和增加气化核心的双重功效，进一步增强对流换热性能并降低沸腾起始温度。

国内外诸多学者对带肋微通道内独特的气液两相流动和相变传热机制展开了研究。Lie等[6]通过实验对比了制冷剂FC-72在带肋微通道内两相流和单相流的换热效果，发现两相流动的换热系数远大于单相流。Krishnamurthy和Peles[7]进行了并排带肋微通道内的流动沸腾实验，在通道内沿流动方向分别观察到泡状流、泡状-弹状混合流、环状流等不同流型，并且发现带肋微通道的流动沸腾换热效率明显高于光滑微通道。Guo 等[8]分别研究了光滑表面和柱状微结构换热面的流动-喷射复合式强化沸腾换热。结果表明得益于柱间的微对流运动，相比光滑表面，相同工况下柱状微结构的临界热通量（CHF）更高，沸腾稳定性得到提升。Qu 等[10]研究了进口流体过冷度对微通道内流动沸腾的影响，发现当过冷度增大，低含气率区域的壁面传热系数会提高，但是在高含气率区域传热变化不明显。为了深入理解流动沸腾换热的诸多物理机制，学者们[11-19]采用数值模拟的方法对该过程进行研究。杨朝初等[11]模拟了微气泡的生长过程，发现气泡一旦成核后会迅速吸热膨胀，并在气泡与通道壁面之间形成薄液膜。弹状流是微通道内非常重要的两相流型，具有良好的传热传质特性[12]。Magnini等[13]基于Hertz-Knudsen 理论和高度函数法构建了高精度的气液相变模型，对微圆管内弹状流气泡的流动沸腾过程进行模拟。结果表明微通道内两相流的壁面传热系数比相同条件的单相流高出30%，并发现薄液膜的蒸发导热是壁面传热强化的主要原因。徐进良等[14]模拟了高热通量下微圆管内触发种子气泡的沸腾传热，并分析了不同气泡触发频率对壁面冷却的影响。从以上文献可以看出，带肋微通道具有良好的沸腾换热性能，但是对流动沸腾过程中所涉及的气泡运动规律以及肋表面的换热特性还需进一步研究。

因此本文采用数值模拟的方法，构建气-液流动和相变模型，对微通道内气泡绕流加热方肋这一基本的流动沸腾问题进行研究。以气泡行为、流场演变以及方肋各表面的传热特性变化作为主要研究对象，通过分析初始气泡体积及入口雷诺数Re 等不同参数的影响，探究带肋微通道内沸腾气-液两相流动及相变传热传质的规律，进而为实验研究提供指导。

1 模型建立

1.1 物理模型

本文构造一个二维微通道物理模型，如图1所示。微通道长度L=4mm，宽度D=0.2mm，通道壁面绝热。边长B=0.08mm 的方柱设置在通道正中心，方柱壁面为恒定热通量q=5×104W/m2加热。各壁面均为无滑移速度条件。微通道入口通入饱和R113制冷剂（Tsat=323.15K），其标准大气压下各相物性参数见表1。入口设置充分发展速度条件，出口为标准大气压出口条件，并在距离入口0.5mm，靠近通道上壁面的位置初始化一个气泡。

以D为特征长度，定义入口雷诺数Re[式(1)]。

式中，u为入口平均速度；μl为液体动力黏度。定义局部传热系数hx[式(2)]。

图1 几何模型和边界条件

表1 一个标准大气压下R113饱和状态的物性参数

式中，Tw,x为方柱表面局部壁温。

定义局部努赛尔数Nux[式(3)]。

式中，λl为液相导热率。

1.2 网格划分

本文采用ICEM 软件划分网格，部分计算区域网格如图2所示。为了保证计算域内气液两相流动与沸腾传热特性的准确性，对微通道壁面以及加热方柱周围进行了网格加密处理，以确保对薄液膜的刻画。为了验证网格的无关性，在此基础上将原有网格数量提高4倍，方柱换热系数结果与原网格算例相比差异小于1.5%，因此可认为结果是网格无关解。

1.3 控制方程和算法设置

本研究中的数值模型基于如下假设：

（1）气液相变温度恒定为一个标准大气压下的饱和温度；

（2）气相和液相的热物性均设置为该饱和温度下的常值；

（3）不考虑液体黏性加热；

（4）忽略重力作用。

以上假设的根据是：①微通道内压降小于5kPa，其对应的液体相变温度变化小于5%[21]；②最大过热度低于13K，其对应的气相和液相的热物性变化均小于5%[21]；③最大Re 下，黏性耗散对流体的加热仅为壁面加热的0.1%[22]；④通道封闭数Co大于5[Co=(σ/gΔρD2)1/2]，重力作用被抑制[23]。

本文使用的流体体积函数（VOF）模型是一种在固定欧拉网格下追踪气液交界面的方法。气相和液相用网格体积分数a 来表示：a=1 为气相，a=0为液相。气液界面处0＜a＜1，物性为各相物性体积分数的加权平均值。气相及液相的质量方程如式(4)、式(5)所示，气液两相共用的动量、能量守恒方程如式(6)、式(7)所示。

式(6)中的Fs表示作用于气液界面的表面张力，通过连续界面力（continuum surface force，CSF）模型将表面力转化为体积力[式(8)]。

式(7)中的能量源项Sh基于Pan 等[18]提出的“饱和界面体积”相变模型，该模型中认为蒸发过程中气液界面处维持在饱和温度，因此计算可以不依赖于经验系数，其表达式为式(9)。

对应的式(4)、式(5)的质量源项表达式为式(10)。

本文使用基于有限体积法的Ansys Fluent 18.0离散型非稳态求解器，流动模型为层流模型。流场采用PISO（pressure implicit splitting of operator）算法求解耦合的压力-速度方程，对时间项采用隐式格式离散。标量梯度的求解采用“基于节点的格林-高斯”法（Green-Gauss node-based method），采用VOF方法中的“结构重构”方法（geo-reconstruct）追踪气液界面。压力插值采用PRESTO 算法（pressure staggering option），动量和能量方程的对流项均采用三阶迎风格式进行离散。为了避免数值振荡，通过设置CFL 条件实现对时间步长的控制，Courant 数设置为0.05，模拟过程中最大的时间步长为10-7s。质量、动量和能量方程的残差分别控制在10-6、10-8和10-10以下。

图2 计算网格划分

1.4 模型可靠性验证

为验证计算模型的精确性和可靠性，本文首先用R113工质对经典的一维Stefan相变问题[24]进行模拟。图3 给出了气液界面位置随时间的变化曲线，可以看出，本文的模拟结果与解析解吻合良好。

图3 相界面位置随时间变化

其次，Han和Shikazono[25]实验研究了微通道内弹状流的薄液膜厚度变化，本文针对液膜厚度与文献结果进行对比。在Ca=μu/σ=0.0029 时，测量值为10.2μm，模拟值为10.7μm；Ca=0.044 时，测量值为63.2μm，模拟值为66μm，模拟值与实验值比较接近。

最后，本文模拟了圆管微通道内弹状气泡的流动沸腾，并与文献[13]中基于Hertz-Knudsen 理论的蒸发模型进行对比，结果如图4、图5 所示。可以看出本文的气泡位置和局部传热系数均与文献结果吻合较好，说明本文采用的模型和方法能够准确模拟微通道内的相变流动问题。

图4 气泡前端位置随时间变化

图5 壁面传热系数沿轴向分布

2 计算结果及分析

2.1 气泡动态及流场变化

基于上述模型，本文通过改变入口雷诺数（Re=60～360）和气泡初始大小（d=0.04mm、0.08mm、0.12mm）进行了一系列模拟，以研究流体流速和气泡体积对方肋微通道内两相流动和强化传热的影响。

图6 显示的是初始直径d=0.08mm 的气泡在不同Re 下，流过加热方柱时微通道内速度场、温度场以及气泡形状的演变过程。其中，气泡的轮廓在速度云图中用黑线标识，在温度云图中用白线标识，带箭头的线为流线。

从图6左侧图片可以看到随着Re增大，受流体惯性力作用，气泡变形程度加剧。同时分析速度场和流线图，Re 为60 和120 时方柱后表面的尾涡区内有两个对称的涡旋；而当Re提高到240后，尾涡区内的涡旋开始周期性地脱落，增强了下游流体的扰动。当气泡从方柱与壁面的狭缝穿过时，气液界面扫掠方柱上表面的过热液体，形成一层薄液膜，挤压热边界层。对比气泡进出方柱前后，尾涡区的流场和温度场发生了明显改变。

图7 分别是四种Re 下，气泡的当量直径deq随气泡质心轴向位置xc变化的曲线。其中，气泡的当量直径deq定义为同等体积下球形气泡的直径；气泡质心的轴向位置xc可用式(11)计算。

式中，N为计算域的网格数，xi是第i个网格的轴向坐标，ai和Vi分别为该网格的气相体积分数与网格体积。

图6 气泡在4种不同Re下，通过加热方柱过程中三个位置的绝对速度及温度云图

图7 不同Re下气泡当量直径deq随质心xc位置的变化

从图7可以看出，当气泡到达方柱中心时（xc=2mm），气泡体积开始陡增，说明此时相变蒸发已经开始。从曲线斜率的变化可以看出，受过热液体分布的影响，气泡在方柱和尾涡区域（xc=2～2.5mm） 体积迅速增长，而到通道下游（xc＞2.5mm）体积增长则不明显。此外，气泡在低Re条件下的增长速度要远远大于高Re 的情况，尤其是在尾涡区内。出现该现象的原因主要有两点：一是高Re 下，气泡流速快，在过热区停留时间短，蒸发吸热过程不充分；二是在高Re 下，流场以涡旋脱落的方式将尾涡区过热液体带向通道下游，对冷热液体进行混合，降低了尾涡区液体的过热度，从而抑制了该区域的相变蒸发。

2.2 换热效果分析

2.2.1 入口雷诺数Re对方柱传热的影响

图8对比了4种Re下，方柱四个加热面平均传热强化率随时间的变化。其中，平均传热强化率定义为(htp-hsp)/hsp，反映了两相流下表面平均传热系数htp与相同条件下单相流hsp的相对大小。定义量纲为1 时间τ=t/(B/u)，当气泡前端靠近方柱时τ 设定为0 时刻。从图中可以看出，当气泡穿过方柱时，方柱4个表面的传热均受到不同程度影响。其中上表面的传热强化最为显著：在Re=60下，最大的换热系数超过单相流的4.5 倍。这是因为当气泡穿过方柱时，气泡在表面张力的作用下挤压周围液体，在方柱的表面和气泡界面之间形成一层薄液膜（图4）。此时薄液膜内部流动缓慢，薄液膜的蒸发导热成为了主要换热方式，其热阻较小，强化传热显著[15]。而随着Re 增大，流体惯性力开始占据主导，液膜厚度变厚，导热热阻增大，导致强化传热效果迅速下降。方柱前表面和后表面也形成部分液膜，因此具有类似变化规律。而对于方柱下表面，虽然未与气泡直接形成薄液膜，但由于气泡堵塞上通道而使下通道的液体流动加快，增强了下表面的对流传热强度。

图8 方柱各表面两相流传热强化率随量纲为1时间τ的变化

进一步对方柱4个表面的传热强化率做积分平均，图9(a)比较了4种Re下方柱整体传热强化率随时间的变化；对气泡前端靠近方柱直到气泡尾端离开方柱这段时间进行积分平均，图9(b)统计了方柱的时均-整体传热强化率与Re的关系。从这两个图都可以看出，在气泡流经方柱的过程中，Re越小，气泡流动沸腾对方柱整体传热强化效果越明显。

图9 气泡流经方柱过程方柱整体换热强化率的变化

2.2.2 气泡初始直径对方柱传热的影响

图10 比较了Re=120 下3 种量纲为1 初始直径的气泡（= deq/D），其质心位置xc恰好位于方柱中心位置（即xc=4mm）时气泡的形状、方柱壁面过热度以及对应的局部努赛尔数Nux的分布情况。

图10 不同初始体积的气泡穿过方柱的对比

3 结论

本文基于VOF 方法，采用“饱和界面”相变模型，以R113 为工质对带肋微通道内气泡流动沸腾进行了二维数值模拟，分析了气泡体积和入口雷诺数Re 对气泡增长速率、液膜厚度和方柱表面传热等重要传热传质参数的影响规律。得到如下结论。

（1）气泡穿过方柱时，在气泡界面和方柱表面之间形成一层薄液膜，挤压温度边界层；气泡对尾涡区的流动结构存在扰动作用。

（2）气泡的蒸发相变主要发生在方柱和尾涡区内，且Re越小，气泡的体积增长越快。

（3）薄液膜的蒸发导热是传热强化的主导机制，气泡的流动沸腾对方柱表面换热有明显的强化作用；对于相同体积的气泡，随着Re 增加，液膜厚度逐渐变厚，强化传热作用迅速减弱。

（4）对于相同Re，气泡体积越大，所形成的薄液膜区覆盖面积就越大，液膜蒸发带来的强化换热也越显著；而小气泡的两相流换热效果同单相流相比并无明显差异。

符号说明

B—— 方柱边长，mm

Co—— 通道封闭数，(σ/gΔρD2)1/2

Ca—— 毛细数，μu/σ

cp—— 比定压热容，J/(kg·K)

D—— 通道宽度，mm

d—— 气泡直径，mm

Fs—— 体积表面张力，N/m3

g—— 重力加速度，m/s2

hfg—— 气化潜热，kJ/kg

hx—— 对流换热系数，W/(m2·K)

L—— 通道长度，mm

Nux—— 局部努赛尔数

n—— 单位法向量

p—— 压强，Pa

q—— 热通量，W/m²

Re—— 入口雷诺数，Re=ρluD/μl

Sm—— 质量源项，kg/(m3·s)

Sh—— 能量源项，W/m3

t—— 时间，s

U—— 速度矢量，m/s

u—— 轴向速度，m/s

a—— 体积分数

λ—— 热导率，W/(m·K)

μ—— 动力黏度，Pa·s

ρ—— 密度，kg/m3

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

τ—— 量纲为1时间，tu/B

x—— 轴向坐标，mm

下角标

c—— 质心

eq—— 当量

l—— 液相

sat—— 饱和状态

sp—— 单相流

tp—— 两相流

w—— 壁面

v—— 气相