脉动热管不同热负荷下管内流型及传热特性分析

黄晓明，卢晓剑，许国良，陈新涛2，段 洋

(1.华中科技大学能源与动力工程学院,武汉 430074；2.张家港科瑞电气有限公司,苏州 215624)

0 引 言

脉动热管(Pulsating Heat Pipe,PHP)自20世纪90年代初被提出以来，就备受国内外研究者的关注。相比于传统的热管，脉动热管不需要毛细吸液芯，利用液塞在蒸发器和冷凝器之间的自激振荡来实现高速热传递，具有结构简单、体积小、热传输能力高等优点，被认为是应用于微小空间、高热流密度条件下的一种极具潜力的高效传热元件[1]。

尽管结构十分简单, 脉动热管的运行机理却非常复杂，是集显热传热、潜热传热、膨胀功于一体的涉及多种机制、多个参数的汽-液两相流系统。大量的实验研究表明，脉动热管的传热特性受多种因素影响，如管径、充液率、工质热物性及工作条件等[2-3]。

苏磊等[1,4]将脉动热管简化为一个热力系统，分析了脉动热管的热力学平衡问题，揭示了能量与系统稳定性的关系；汪健生等[5]研究了脉动热管内的气液两相流动及工质所涉及的沸腾情况，并通过数值模拟，分析了固定热负荷下加热段与冷凝段长度比值对脉动热管启动和运行过程的影响；邵帅等[6]将脉动热管简化为一维直管，通过建立脉动热管多气塞液塞模型对以氮为工质的脉动热管进行了研究，分析了管内振荡和传热特性，并得出了脉动热管内工质在温差驱动下的周期运动。不过其采用的多气塞-液塞模型采用了较多的简化假设，不能反映复杂的两相流变化对系统传热性能的影响。Sakulchangsatjatai等[7]依据可视化试验的结果，观测到了气液两相之间的转化过程，并基于连续介质力学分析方法对管内两相流过程进行了数值模拟。Hanlon等[8]研究了蒸发段与冷凝段温差、气液塞初始温度、弯管数量等因素对脉动热管内部工质振荡的影响，揭示了脉动热管内部的流体动力学情况，但并未进一步探讨动力学特性与系统传热性能之间的定量关系；Ma等[9]对管内气液两相流的运行机理进行了分析，指出工质振荡行为对脉动热管传热性能有很大影响，但在二者的定量关系方面也未做进一步探讨。

目前，人们对液塞在脉动热管中自激振荡运动规律的研究还不够完善，较多的简化条件使得模型与系统实际传热特性之间难以建立直接联系。且大多数基于两气塞-液塞振荡模型的研究都只停留在一个传热单元的比较分析。根据可视化实验和模拟结果分析[10]，不同热负荷条件下，脉动热管内流体流动会发生弹状流到环状流的流型变化，这种流态的变化对蒸发段管壁与工质之间的传热系数有很大影响，也会对热管整体流动和传热特性产生很大差异。因此，在建立脉动热管机理模型时必需考虑流体流型的变化及其带来的传热特性差异。

综上所述，文中将复杂的气液两相传热传质过程融入到动量和能量平衡机制中，对传统的两气塞-液塞振荡理论模型进行了改进。在新模型的基础上，探讨了热负荷、管径和传热单元数等因素对脉动热管流动和传热特性的影响。

1 物理模型和数学模型

脉动热管的结构如图1(a)所示。热管分为蒸发段、绝热段和冷凝段，汽塞与管壁之间有液膜层，液膜厚度与距液塞的远近有关。在绝热段内，工质无吸热和放热过程，建模时可忽略绝热段的影响。文中通过取脉动热管的一个U型单元，构建的物理模型是一个垂直放置的U型圆管微通道，其两端密封(如图1(b))。将微型通道视为直管(如图1(c))所示，考虑流动过程中重力的影响，忽略工质在绝热段发生的物性变化，忽略弯管处的压力损失[11-12]。位于管道两端的每个蒸发段长度为Le，蒸发段管壁温度保持在Te。冷凝段部分位于两个蒸发段之间，长度为Lc，冷凝段管壁温度保持在Tc。液塞长度为Lp，与脉动热管的工质填充率有关。用xp表示管内液塞的运动位移。当液塞向右移动时，xp为正；液塞向左移动时，xp为负；液塞最左端与冷凝段最左端重合时，xp为0。

图1 物理模型

当液塞的初始位移xp为正时，左侧汽塞部分进入冷凝段，发生冷凝，汽塞压力Pv1减少，热管右侧的液膜蒸发导致右侧汽塞压力Pv2增大，左右汽塞的压力差导致液塞向左移动。当xp变为零时，两个蒸汽塞中并没有蒸发、冷凝，但由于惯性，液塞继续向左运动。当xp变为负时，左侧汽塞压力增大，右侧汽塞压力减小，压差方向改变，液塞运动到最左侧后开始向右运动。两个蒸汽塞交替蒸发和冷凝，压差方向交替变化，维持液塞的振荡运动，实现热量的传送。

模型中做出以下假设：

(1)设各结构单元之间的传热传质过程相似；

(2)个蒸发段的汽塞可视为理想气体，液塞是不可压缩流体；

(3)汽塞和液塞之间面上发生质量和能量交换全部是由于相变引起的；

(4)壁上的剪切应力与流体的流动状态有关。液塞运动方向改变时，剪切应力方向改变。

1.1 液塞运动方程

对液塞进行受力分析：运动过程受到两侧汽塞压差产生的推力、液塞两侧高度差产生的重力、液塞与管壁之间的剪切力。根据牛顿第二定律，当液塞由左向右移动时，有公式(1)：

(1)

式中，A为脉动热管的横截面积；τp为液塞与管壁的剪切应力，可用公式(2)表示：

τp=0.5Cρu2

(2)

式中,u为液塞的运动速度，C成为粘滞系数，与液塞的流动状态有关。根据液塞的运动速度计算出雷诺数Re(公式(3))后，可以由公式(4)计算得粘滞阻力系数C的值。

(3)

其中,ν为液塞的运动粘度系数。

(4)

综上，将公式(1)简化为：

(5)

二阶微分方程，可以通过四阶龙格库塔法求解得到对应时刻的位移值xp。

1.2 汽塞能量方程

蒸汽塞的能量的变化量，其值为液膜蒸发进入汽塞的能量与蒸汽塞对外界做工能量的差值。左右侧汽塞的能量方程分别为公式(6)和公式(7)：

(6)

(7)

式中，hfg为该时刻对应气塞温度下的汽化潜热值，cv为定容比热容值，Tv1、Tv2分别做左、右侧汽塞温度，mv1、mv2分别为左右侧汽塞质量。

假设两个蒸发段的汽塞遵循理想气体定律，写出蒸汽的状态方程(8)和(9)：

(8)

(9)

方程(8)和(9)分别对时间求导得到方程(10)和(11)：

(10)

(11)

将方程(10)代入方程(6)，方程(11)代入方程(7)，将得到的微分方程进行积分，在给出汽塞初始压力Pv10、初始位移xp0、初始质量mv10的情况下，可以得到气塞压力、位移和质量的关联方程(12)和(13)：

(12)

(13)

1.3 汽液塞质量交换

蒸汽塞的质量变化与液塞所处的位置有关，假设汽塞与液塞的质量交换只发生在汽液交界面上。蒸汽塞质量变化可以由公式(14)和(15)描述[12]：

(14)

(15)

1.4 两相沸腾换热分析

已有研究表明，不同传热负荷条件下，脉动热管内工质会呈现出不同的流型。当其处于稳定运行时，主要流型为弹状流与环状流。文中借助Chen[13]提出的经验公式来计算he，将脉动热管内传热视为沸腾换热和强制对流换热两种方式的叠加，即：

he=EhFc+ShNB

(16)

其中，hFC是强制对流传热系数，hNB是泡核沸腾传热系数。E为对流强化系数，S为泡核沸腾抑制因子。

强制对流传热系数与流体热传导能力、流动速度及含气率有关，文献[14]推荐经验公式：

(17)

式中，k1、μ1和cp1分别为液体的导热系数、动力粘度和定压比热容；Pr1和Re1分别为以液体物性为基础计算的雷诺数和普朗特数；d为脉动热管管内径，G为流体质量流流速，kg/(m2·s)；F两相流动参量，反映两相流动对传热的影响；x是质量流含气率，为气相的质量流量Mg和气液混合物的总质量M之比，即：

(18)

其中，Ag是气相流通界面积，wg是气相流速。

式(17)中两相流动参量F由Chen经过实验确认，F与Xtt的关系为：

(19)

Xtt为气液两相均为湍流时的Martinelli参数，根据质量流含气率x确定：

(20)

hNB采用Liu-Winternton 改进泡核沸腾关系式[15]：

(21)

其中，M为流体分子量，水为18；Pr=p/pcr为折算压力，q为热流密度W/m2。

泡核沸腾因子S表达式：

(22)

(23)

对于弹状流型下对流强化传热系数E，采取Dittus-Boelter关系式[16]计算得到：

E=1+1.37(1/Xtt)0.86

(24)

环状流下对流强化传热系数E，需要考虑沸腾剧烈程度，增加沸腾项：

E=1+24 000Bo1.16+1.37(1/Xtt)0.86

(25)

其中，Bo为沸腾数：

(26)

式中，il-g工质的汽化潜热，(J/kg)。

2 求解过程

基于MATLAB平台对上述方程采用数值求解，其主要求解过程如下：

(1)给定蒸发段和冷凝段管壁温度、环境温度、热负荷等边界条件，对管内液塞的速度、位移及管内气塞的物性参数变化进行预判断，并根据预判断确定管内两相流型，并根据流型计算蒸发段壁面传热系数；

(2)据脉动热管内液塞运动方程、气塞运动方程和气液塞质量交换关联式，求解两侧汽塞的初始质量；

(3)利用四阶龙格库塔法[17]，求解微分方程(5)，得到液塞的位移xp和速度u；使用牛顿-拉夫森算法(Newton-Raphson)求解方程(12)和方程(13)中的左右侧压强校验值Pv1和Pv2。

(4)蒸发换热量、流型参数进行校验。满足假设后认为过程收敛，并输出系统热分析所需要的参数，进行下一步的系统换热分析。

3 结果与讨论

针对图1的5单元闭合回路脉动热管，设定如表1所示基本参数，进行模型验证和分析。

表1 脉动热管基本物理参数

续表1

物理参数数 值冷凝端长度Lc0.4m初始位移xp0热负荷Q15, 25, 40, 55, 70, 80W环境温度25℃管内工质水充液率50%传热单元数N4,5,6

3.1 模型验证

为验证模型的准确性，文中将数值计算结果与文献[18]的实验测量数据进行了比较。文献实验研究所采用的脉动热管内径为2 mm，传热单元数5，其它参数与表1中数据完全一致。二者的比较如图2所示。由图可以看出，模型能够很好地预测脉动热管传热特性的变化，在趋势和数值上都与实验数据非常吻合。基于图2的理论计算结果，我们还可以预测到：当脉动热管处于低负荷运行，管内流型主要以弹状流为主时，热管热阻较大；随着热负荷增大，脉动热管内的流型逐渐演变为环状流，此时热管热阻较小，且随着热负荷的增大，热阻会逐渐降低，但变化越来越小。

图2 PHP模型计算热阻与试验热阻比较(N=5,d=2 mm)

3.2 热负荷影响

脉动热管运行热负荷对脉动热管管内工质流型和工质与管壁的总传热系数具有较大的影响。不同管径下，热负荷对热管传热特性的影响也有不同的表现。由于管径改变，为保证相同的管内工质流通体积，需要相应改变加热段和冷凝段的长度。比较1.5、2、2.5 mm管径下，运行热负荷与脉动热管管内流型和传热系数he的关系，结果如图3所示。

图3 不同管径下，热负荷与传热特性关系(N=5)

对脉动热管蒸发段内两相流流动和传热特性随热负荷的变化情况进行了分析，结果如图3所示。

图3(a)反映了蒸发段内两相流体的分液相动压头ρljl和分气相动压头ρgjg随热负荷的变化。两个纵坐标的计算公式如下：

(27)

依据文献[19]给出的两相流型判定图，认为当ρljl≤104，且ρgjg≤103时，管内以弹状流为主；当ρljl≤104，且ρgjg>103时，管内以环状流为主。由图3(a)可以看出，流型在较低的热负荷下易形成弹状流，随着负荷的增高逐渐演变为环状流。在相同热负荷下，不同的管径流型也会有不同的表现，如同为25 W热负荷，1.5 mm管径内以环状流为主，而2 mm和2.5 mm管径则以弹状流为主。这是因为小管径脉动热管蒸发段工质汽化产生的气泡更容易聚集，气泡成长过程排开液态工质体积并不管合并，容易产生环状流。

图3(b)可以看出，当脉动热管管内从弹状流演变至环状流时，对流换热系数he会有一个比较大的陡增，且随着热负荷的增大， he继续增大。2 mm管径时，40 W热负荷下环状流he是25 W热负荷下弹状流he的近5倍。各管径不同流型下的he差异非常大，且都随着热负荷在逐渐变化。

结合图3(a)和图3(b)可知，两相流流动形态与热负荷和管径等多种因素有关。不同流态对应的传热特性不同，有必要在模型中加入相关的考虑。

3.3 PHP振荡特性分析时

针对管径为2 mm的脉动热管，不同热负荷下脉动热管液塞运动情况如图4所示。

由图4(a)可以看出，低负荷条件下(25 W)，脉动热管内以弹状流为主时，液塞振荡幅度较大，比环状流最大振幅大0.08 m，频率也较低；而较高负荷条件下(55 W)，管内以环状流为主时，管内的液塞振荡幅度较小，频率却高达10 Hz，比弹状流的振荡频率高出2 Hz左右，振荡更为剧烈。随着热负荷的进一步增长(80 W)，振荡幅度也进一步减小，振荡频率进一步增大，但因为没有流型的差异，振幅和频率相比55 W的情况并没有太大变化。

图4 PHP振荡特性图(N=5,d=2 mm)

图4(b)和图4(c)是左侧汽塞温度和压力随热负荷的变化，可以看出，25 W时温度和压力变化幅度较大，这是由于低负荷下的弹状流振荡程度不够剧烈，引起的热量在蒸发端聚集，未能及时传递到冷凝段，引起的加热段温度和压力的升高。对于较高负荷，管内流型演变为弹状流后，频率升高，振荡剧烈，温度和压力的变化速度加快，变化幅度降低。图4(c)还显示，热负荷越大，气塞的压力被整体抬高。

3.4 几何参数对PHP传热特性影响

图5 直径d、传热单元数N对PHP整体热阻影响(控制热管工质的流通总体积一定)

PHP的传热特性受到许多物理参数的制约，如管径d、传热单元数N等。为探究管径对脉动热管整体运行热阻的影响，针对图1所示5弯头PHP，选取1.5、2、2.5 mm管径的PHP传热运行情况，比较各管径下热管的整体热阻值，如图5(a)所示。

由图5(a)可以看出，脉动热管内工质流型以弹状流为主时，热阻较大，而演变为环状流后，热阻会有较明显的降低。所以针对低负荷运行的工况，可以选择管径较小的脉动热管。随着负荷的增加，当各管径内工质运动都演变至环状流为主时，三者的整体热阻差异不明显，管径打的热阻会略低一些。针对更高负荷的工况，不宜选择小管径的脉动热管。

为探究不同传热单元数N对脉动热管弯头数目对运行热阻的影响，选取N为4、5、6，管径d都为2 mm的PHP传热运行情况，比较各传热单元数N下热管的整体热阻值。需要控制加热和冷凝段的总长度这一变量恒定，故在改变N时，需要相应改变加热段和冷凝段的长度。结果如图5(b)所示。

由图5(b)可以看出，在相同热负荷下，随着传热单元N的增加，脉动热管的整体热阻会减小，传热特性会提升，但提升效果不明显。特别是处于低热负荷区域，三根不同数目的脉动热管传热特性差异较小。由图中也可以看出，N的数目变化对流型变化趋势的影响并不大。

4 结束语

文中对脉动热管两气塞-液塞振荡理论模型进行了完善，将复杂的两相流传热传质过程与脉动热管内动量和能量平衡机制相结合，使模型能够更好地反映两相流动和传热特性对脉动热管整体热传输的关系。根据文中的研究结果，得到如下的结论：

(1)通过与已有文献的实验相比较验证了模型的准确性。计算结果与实验数据的吻合，不仅反映在趋势上，还反映在数值上。同时，文中的模型还可以对流态进行预测和传热特性进行分析。

(2)低负荷下，脉动热管内流体以弹状流为主，传热系数较小，振荡运动为低频高幅；随着负荷的增大，流态将从弹状流转变到环形流，其转变的临界热负荷与内径有关。环形流时，传热系数有明显增大，振荡运动特性为高频低幅。

(3)从影响因素看，流态与热负荷的关系主要取决于内径，与弯头数目关系不大。但弯头数目的增大，在同样条件下，会使系统热阻略有下降。