单环路液氢温区脉动热管高充液率工况计算流体动力学（CFD）模拟

徐金柱，焦波，孙潇，王芳，甘智华

（1 哈尔滨理工大学荣成校区机械工程系，山东荣成264316；2 浙江大学制冷与低温研究所，浙江杭州310027）

随着电子芯片的小型化，其局部热通量已经高达102～103W/cm2[1]，严重威胁到设备的安全运行。为了解决这一问题，亟需开发一种高效的传热元件。脉动热管（pulsating heat pipe，PHP）具有体积小、结构简单、传热效率高等优点，一经提出就引起了研究者们的注意。它由一根毛细管弯折而成，充注完成后由于毛细力的作用工质在管内形成随机分布的气液塞，工作时，气塞会由于蒸发和冷凝的作用而膨胀或缩小，推动工质在管内形成振荡或循环流动。与传统热管相比，PHP启动更快，传热温差更小，在低温超导[2]、空间应用[3]、余热回收[4]等领域具有广阔的发展潜力。但PHP内部为复杂的气液两相流，目前仍无法完全理解其运行机理[5]，并且影响PHP传热性能的因素较多，增加了其研究难度。

工质对PHP 的传热性能有很大影响，它不仅决定了PHP 的管道临界直径，而且不同工质PHP的换热能力有很大不同[6]，工质的选择取决于PHP的运行温度和应用范围。徐冬等[7]总结了近些年来低温脉动热管的实验及理论研究进展，他们指出低温温区PHP 的研究起步较晚，影响传热性能的参数尚未被完全确定。随着低温超导材料（临界温度39K）的发现和应用，液氢温区PHP的高效传热引起了研究人员的关注。2011 年，Natsume等[9]对氮、氖、氢低温PHP 的传热性能进行了测试，结果表明氢工质PHP的有效热导率高达500～3500W/(m·K)。2014 年，甘智华等[10]设计了针对液氢温区PHP的实验台，以探索多种参数（弯头数、加热功率、充液率、倾角等）对PHP 传热性能的影响。他们[11]测试了充液率34.2%情况下脉动热管的运行过程，随后，变充液率和变加热功率下PHP的性能得到了实验研究[12-13]。2019年，孙潇和甘智华等[14]在实验中将脉动热管的冷凝段温度分别维持在不同数值，分析了不同的冷却温度对液氢温区PHP传热性能的影响，发现冷凝段温度的提高可以增加PHP 的传热效率，并不会对PHP 的启动造成影响。

另一方面，研究PHP 传热机理的手段也得到了极大发展。随着计算机内存的扩大，计算速度的加快及算法的优化， 计算流体动力学（computational fluid dynamics，CFD）数值模拟被大量运用于分析PHP的流动及传热特性。刘向东等[15]建立了3D 数值模型分析了板式脉动热管内工质的流动方式，他们指出工质在板式PHP 存在三种流动方式：离散的气泡流，气塞流和长气塞流。Pouryoussefi和Zhang[16]建立了一个2D模型对PHP进行了分析，探究了以水作为工质的脉动热管的流动情况。他们指出，对于闭式PHP，当蒸发段的热通量较高时，从蒸发段向冷凝段运动的液塞具有较大的速度，这个速度足够大以至于液塞可以穿越冷凝段到达下一个蒸发段。此时，PHP内的流动即从振荡流变为了循环流。Wang等[17]建立了2D模型对充液率分别为30%、50%和70%的PHP 进行了数值模拟，他们得出，在底部加热模式下，随着加热功率的增加，不同充液率的PHP传热热阻均会减小，而充液率较低的PHP具有更小的热阻。随后他们[18]对单环路水工质PHP 的传热性能进行了CFD 模拟,利用2D 模型，分析了以水为工质的PHP 在加热功率变化时热阻的大小，并将模拟结果与Saha等[19]所做的实验相对比。他们指出，模拟与实验得到的热阻变化趋势一致，并且它们之间的误差为10.01%。在低温PHP 的CFD 模拟中，陈曦等[20]在对氮工质PHP 进行模拟时发现，在PHP 稳定运行后，管内工质的循环流动并不是一直在进行，而是在运行一段时间后会经历一段时间的停滞，然后重新运行。

综上所述，液氢温区脉动热管的实验研究已经积累了一定的研究成果，CFD数值模拟也被广泛地应用到PHP 理论分析中，同时有效性得到验证，但尚未发现针对液氢温区PHP 的CFD 模拟研究，启动阶段及运行阶段脉动热管管内的压力变化，以及压力变化与脉动热管流动方式的关系，在以往的CFD 数值模拟中均未发现相关报道。由于氢的密度、汽化潜热远小于常温工质，饱和压力对温度的导数dpsat/dT高于常温工质一个数量级[14]，导致管内压力随温度的变化更加剧烈，在传热性能及运行机理也与常温脉动热管有所差异，因此开展低温液氢温区脉动热管的研究对于全面理解脉动热管传热机理与流动规律具有重大意义。由于低温工况下实验难度大、耗时长，可视化实验尤其困难，利用CFD模拟可以获得PHP 的传热性能和内部流动特性。基于此，在本文作者课题组已有的研究基础之上[5,10-13,21]，本文针对单环路液氢温区PHP建立了二维数值模型，对80%充液率工况进行了模拟，并与实验对比验证了模型的有效性，分析了在启动过程中PHP 内工质的流动与管内的压力、温度的关系，探索了稳定运行阶段PHP 中温度及压力的振荡规律，并讨论了蒸发段的加热量对PHP 传热性能的影响，旨在发展液氢温区PHP理论分析模型，为指导工程应用提供依据。

1 控制方程

在本研究中，采用VOF 方法追踪PHP 中气液界面。VOF方法[22]适用于两种或两种以上不相溶的流体。因此，根据计算单元中体积分数的大小可以将控制体分为两种情况：一种是纯流体，αk=1；另一种是混合流体，αk≠1。在PHP 中只存在气液两相，当控制体中只存在气相或液相时，控制方程与单相流的控制方程相同。当控制体中存在相界面时，控制体中气相与液相的体积分数之和为1，如式(1)所示。

流动及传热的控制方程见式(2)～式(5)[23]。

连续性方程

动量守恒方程

能量守恒方程

动量方程和能量方程所求解出的速度场和温度场被气相和液相流体共享[24]，因此混合相的物性被用于动量方程与能量方程中，见式(6)。

在PHP 运行过程中，发生在气液界面上的传热和传质是引起界面变化根本原因[25]。为了计算相变过程中传质量和传热量，Lee 等[26]提出了一个传质模型，开发了一套完整的CFD 代码，并将其连接到控制方程中。确定了传质量后，相变过程中的传热量由传质量与汽化潜热的乘积而计算[式(7)～式(10)]。

2 物理模型和边界条件

PHP管径必须足够小，才能使工质在管道中呈气液塞分布，此时，脉动热管的临界直径Dcrit需满足以下条件[式(11)][27]。

密度和表面张力是关于温度的函数，那么PHP的临界直径同样依赖于工质的运行温度，图1表示的是氢工质脉动热管的临界直径与运行温度的关系，在实验中取脉动热管的内径为2.3mm，此时对应的最高工作温度为28.6K。当运行温度低于此温度时，该值均小于临界直径，满足PHP运行条件。

图1 临界直径与运行温度的关系

由式(11)可知，重力和表面张力在内径的选择上起到了重要作用，由于本文所建立的是二维模型，相当于将物理模型沿纸面拉伸1m[24]。为了比较两者的差别，如表1所示，对二维模型与三维模型的弯液面进行了分析，假设微元体高度为ΔL。

从以上量纲分析可知，在相同表面张力的作用下，二维模型的内径为三维模型内径的一半。因此在模拟中内径为实验内径的二分之一至1.15mm，汪健生等[28]在对水工质脉动热管模拟时采用了相同的内径处理方式。

表1 二维模型与三维模型的区别

由以上分析，本文建立了单环路PHP的二维数值模型，如图2(a)所示，垂直放置并采用底部加热模式。蒸发段、绝热段和冷凝段的长度分别为51.5 mm、100mm和51.5mm。工质为氢，充液率为80%。

图2 物理模型

在本研究中，作出了以下假设：①忽略壁面厚度；②气体密度遵循理想气体状态方程，其余物性为温度的函数；③饱和温度为压力的函数；④流动为层流。为了保证网格的质量和改善计算精度，在模拟中采用ICEM 生成结构化网格。图2(b)展示了部分蒸发段的网格示意图，近壁面处的网格被加密。

根据实际的物理过程，整个模拟过程分为3个阶段。

（1）初始阶段 实际过程中工质充注到PHP后，由于毛细作用，工质在PHP 内会随机形成交替分布的气塞和液塞。在模拟时，为了得到这种气液塞随机分布状态，初始化后，将所有壁面边界条件设置为冷凝段温度（19K）。初始时，管内工质的温度也与此相同，因此在初始化时，将管内压力设置为初始温度对应状态下的饱和压力。

（2）启动阶段 管内形成相对稳定的气液塞分布后，蒸发段、绝热段和冷凝段的边界条件分别为：恒热通量、绝热和恒温。维持冷凝段的温度为19K，蒸发段的热通量为加热量与蒸发段的表面积之比。

（3）稳定阶段 当蒸发段温度在小范围周期振荡即认为PHP 已经达到稳定阶段。模拟中在启动成功并且温度达到稳定后，逐渐提高蒸发段的加热量，以分析其对PHP传热性能的影响。

本模拟中，所有控制方程由基于有限体积法的软件Fluent 求解。PISO（pressure implicit splitting of operator）算法用于迭代离散方程，压力的差值格式为PRESTO!格式，动量和能量方程的迭代采用二阶迎风差分格式。模拟中测试结果显示采用默认的松弛因子就可以获得收敛的结果。用于计算的时间步长为10-4s，收敛标准分别是：连续性方程5×10-4，速度10-4，能量方程10-7。

3 结果与讨论

3.1 模型的有效性验证

为了检验网格独立性，验证数值结果的可靠性，采用了4种不同的网格，并计算得出了相应的网格数对应的热阻，如图3所示。考虑数值模拟的计算精度，在本研究中网格数量取为58772。

图3 网格独立性验证（Q=0.5W）

为了验证数值模型的有效性，将模拟结果与文献[14]中的实验结果进行了对比，实验中采用了弯头数N=2的脉动热管，在本文中对单环路脉动热管进行了数值模拟，为了比较两者的传热性能，将加热功率减小到文献中的1/2，以此保证每个弯头的热通量相同。PHP的传热性能用热阻来说明，总热阻的表达如式(12)所示。

图4给出了模拟与实验热阻的对比，可以看出两者具有良好的一致性，最大的热阻误差不超过15%，这是因为在模拟过程中会不可避免地由于几何建模、假设条件等引入数值误差，并且模拟与实验的边界条件不可能完全相同，导致模拟与计算存在一定误差。从图4中可以看出，在充液率为80%的实验工况下，随着加热功率的增加，脉动热管的热阻逐渐减小，但热阻变化率随着加热功率的增加而降低，CFD数值模拟得到了与之相接近的变化趋势。这是因为当加热功率增加后，蒸发段的工质更容易发生核态沸腾，液膜蒸发速度加快，相变速率提高，导致脉动热管内工质的振荡频率提高，从而提高了传热性能降低了热阻。随着加热功率的增加，热阻逐渐趋于稳定。这是由于脉动热管流动需要克服的流动阻力，由于充液率较高，管内工质的流动更接近单相流，热阻随着加热功率的变化更容易达到稳定。

图4 实验热阻与模拟热阻对比

3.2 初始阶段PHP内压力分布特点

当工质充注到PHP 后，由于毛细作用管内的工质会形成随机分布的气塞与液塞，图5为初始阶段管内气体体积分数分布及压力分布云图。从整体图上看，PHP内部的压力存在一个梯度，高压区位于蒸发段，最高压力出现在蒸发段的底部，两个通道沿着冷凝段方向压力逐渐下降，这是由于重力在初始阶段的气液分布与压力分布中起到了重要作用。在局部区域，压力波动下降，高压与低压相互交错，图6展示了轴线方向上[P-Q是图2(a)中蒸发段右侧的直管段部分]压力与对应位置的体积分数的变化。可以发现压力的变化和气液分布存在一定的相关性，高压区域对应着气塞位置，低压区域对应液塞位置。气塞和液塞间存在清晰的相界面，界面之间存在一个明显的压力跳跃[29]，这是由于表面张力和液膜弯液面引起的。初始阶段的气液分布与压力分布为PHP的启动提供了基础。

图5 初始阶段气体体积分数分布云图及压力云图

图6 初始阶段部分蒸发段的压力分布

3.3 启动阶段PHP的流动分析

脉动热管的启动与管内压力的变化息息相关，在启动阶段，蒸发段与冷凝段依旧由于重力的影响而存在压差，但由于液体的蒸发和冷凝使PHP 内压力变化的量级远大于此，因此可以忽略蒸发段和冷凝段的压差，通过监测脉动热管蒸发段内点6#流体的压力[6#的位置如图2(a)所示]，从而得到PHP管内压力随时间的变化曲线，如图7 所示。同时，在模拟中得到了启动阶段管内气体的体积分数变化，可以表明工质在不同时刻的流动状态，由此分析压力变化对工质流动的影响。在本研究中工质经历了两次循环（图7中ABC和CDE）完成启动。在每次循环中，工质都经历了以下3个过程：通道内同时向上流动、单方向循环流动和逆向流动。

图7 PHP内流体压力波动曲线（Q=0.27W）

经过初始阶段的计算，脉动热管内形成了随机分布的气塞和液塞。从图7中A点开始对PHP的蒸发段施加恒热通量边界条件，脉动热管开始进入启动阶段。蒸发段的壁面被加热，内部的流体发生核态沸腾，一些汽化核心随机形成，产生小气泡，原来处于蒸发段的气塞也会随着液膜蒸发而膨胀，由于PHP 蒸发段的加热面是对称布置的，蒸发段两根通道内的气塞会由于膨胀作用同时向冷凝段运动，图8(a)展示了蒸发段右侧通道P-Q 的气塞变化。与此同时，在冷凝段中由于PHP 压力升高导致饱和温度升高，气塞温度低于饱和温度，气塞开始冷凝，逐渐缩小为小气泡，原来的小气泡冷凝为液相，如图8(b)所示。此阶段工质在通道内同时向上流动。

随着蒸发段气塞液膜的不断蒸发，气泡合并，蒸发段中气塞逐渐变长，脉动热管内的压力也不断增加。2.4s时脉动热管内的压力与气体体积分数分布如图9(a)所示，可以发现蒸发段左右两侧通道内的压力明显不同，左侧通道内的压力低于右侧通道的压力。这是因为汽化核心的产生位置与蒸发作用都具有一定的随机性，蒸发段的气塞数量、大小因此会产生差异，气液分布状态发生了变化，PHP蒸发段左右两侧通道出现压差。当压差足够克服重力和剪切力时，管内工质开始向同一方向的流动，如图9(b)所示，在此工况中，PHP内的工质最开始沿顺时针方向流动。而此时PHP 整体的压力仍在增加，直到3.3s附近，气塞即将运动到冷凝段，压力达到最高值（图7中B点）。如图9(b)所示，此阶段工质的流动方式为单方向循环流动。

PHP的压力最高值对应饱和温度的最高值，蒸发段的工质需要达到更高的温度才能发生核态沸腾与蒸发，而此时工质刚刚从冷凝段到达蒸发段，需要较长的加热时间，因此在3.5～4.0s 期间蒸发段几乎全部为液相，无法为工质运动提供必要的压力。由于惯性作用，工质会继续向前流动一段时间后转变流动方向，如图10 所示。随着气塞在运动中冷凝使管内压力降低，饱和温度降低，蒸发段的工质重新开始发生核态沸腾，产生小气泡，压力开始回升（图7 中C 点）。由于工质的逆向流动，使蒸发段右侧通道的工质在蒸发段的加热时间比左侧通道的工质加热时间长，因此它会比左侧通道更早发生核态沸腾，产生更多的气泡，如图10 中4.2s所示。此时，工质的流动方式为逆向流动。

随着加热的进行，蒸发段左右两侧出现核态沸腾现象，同时冷凝段的气塞也相继液化，如图11所示。PHP两根通道中的工质重新开始同时向上流动。由于右侧通道比左侧通道具有更多气塞，工质比上一次更早形成单方向循环流，气塞更快地到达冷凝段，使得管内压力的峰值（图7中D点）低于第一次的峰值（图7 中B 点），并在随后工质改变流动方向，完成第二次循环。随着循环的进行，压力的振荡趋于稳定，PHP进入稳定阶段。

图8 气塞的变化

图9 工质沿顺时针流动开始时压力与气体体积分数云图及体积分数随时间的变化

图10 PHP气体体积分数变化及流动方向转变

3.4 稳定阶段PHP的流动及传热特性

在稳定运行过程中，脉动热管中的流动主要为顺时针循环流，这种流动方式使通道出现上升管和下降管，工质顺时针流动时，左侧通道为上升管而右侧通道为下降管。工质的温度与气体体积分数分布如图12 所示，可以发现上升管中流体的温度始终高于下降管中流体的温度。这是因为下降管中的流体刚经过冷凝段的冷却，流体温度较低；而上升管中的流体在蒸发段经过了充分加热，在上升管发生核态沸腾，产生气泡，气泡相互合并形成气塞，将工质分隔为气液塞分布状态，到达冷凝段后冷凝。

图11 PHP内气体体积分数变化（4.3～5.9s）

工质的温度变化是判断PHP 是否成功启动的重要指标。PHP 内的温度变化如图13 所示，温度振荡的峰值均已在图中标出。

其中，Twe为蒸发段壁面的平均温度；Tev和Tcv分别为蒸发段和冷凝段流体的平均温度，由式(13)和式(14)计算[1#为蒸发段轴线上的点，位置如图2(a)所示]。

图12 稳定阶段PHP内工质的温度分布及体积分数分布

图13 PHP温度随时间的变化（Q=0.27W）

与图7压力变化曲线对比可以发现，蒸发段壁面温度的变化与压力的变化趋势相同，但在启动阶段温度的峰值所对应的时间（2.8s 和5.3s）明显早于压力峰值所对应的时间（3.3s 和5.8s）。开始时，PHP 的温度振荡较大，主要是因为刚开始PHP 中工质运动速度较小，蒸发段的气塞持续蒸发，部分液膜被蒸干，气塞直接接触壁面，导致气塞及蒸发段的壁面温度急剧升高。而当工质开始循环流动后（2.4s后），绝热段的流体进入蒸发段，使蒸发段的壁面温度开始急剧下降，随着工质的不断流动，壁面温度也开始趋于稳定。值得一提的是，蒸发段流体的平均温度则可能会低于冷凝段流体的平均温度。这是由于流体速度较高，从冷凝段流向蒸发段的冷流体在蒸发段尚未被充分加热的原因。

3.5 加热功率对PHP蒸发段温度波动的影响

在模拟中改变蒸发段的热通量，使加热量逐渐从0.27W 增加到1W，得到了PHP 温度随加热量的变化曲线，如图14 所示。在0.27W 加热量下，蒸发段壁面温度有规律地振荡，当加热量升高到0.5W 时，壁面温度小范围升高，随后稳定，而在0.5W 和0.75W 加热功率下，PHP 的壁面温度振幅很小，在功率增加到1W时，振幅突然增加，而当加热量增加到1.25W 后，壁面温度振幅逐渐发散，温差增加，传热逐渐恶化。

图14 脉动热管温度随时间的变化曲线

从图14 中可以发现，在稳定运行阶段，温度的振荡有周期性，对不同功率下稳定运行阶段的温度进行快速傅里叶变换（FFT），得到的信号强度与频率的关系如图15 所示。从图中可以发现各个功率下FFT分析均存在一个显著的峰值，对应的频率如图中所示。这说明PHP 在稳定运行阶段存在显著的周期性，运行频率随着加热功率的升高而升高，而当加热功率达到一定值时，运行频率趋于稳定。通过与热阻的变化（图4）进行对比，可以发现，当加热功率较低时，温度振荡的周期长，频率小，热阻大，随着加热功率的升高，温度振荡频率加快并趋于稳定，PHP的传热热阻随着振荡频率的加快而减小，并最终趋于稳定。

4 结论

本文对单环路液氢温区PHP 建立了2D 模型，利用CFD 数值模拟研究方法，分析了其在80%充液率下的流动及传热特性，研究了不同阶段PHP内压力分布与温度变化。结果表明，CFD数值模拟结果与实验结果误差小于15%，脉动热管的热阻随着加热功率的增加先减小而后不变。初始阶段PHP内气塞与液塞交替形成，由于重力影响，PHP内压力存在明显的梯度，局部位置由于气液黏性的不同，导致压力呈高低交错分布。当有热量输入后，工质经历了两次通道内同时向上流动、单向循环流和逆向流动的过程完成了启动进入稳定阶段，此时PHP的主要流动形式是顺时针循环流，蒸发段壁面温度的变化与压力的变化趋势相同，在稳定阶段呈有规律的周期波动。不同加热功率下，温度振荡的频率不同，它随着加热功率的增加先增大而后趋于稳定。

图15 不同加热功率下稳定运行阶段蒸发段壁面温度的FFT分析（FR=80%）

符号说明

d—— 直径，m

E—— 比焓，kJ/kg

F—— 表面张力引起的体积力源项，N/m2

Fg—— 重力，N

Fs—— 表面张力，N

g—— 当地的重力加速度，m/s2

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

p—— 压力，Pa

Q—— 输入功率，W

R—— 热阻，K/W

S—— 源项，kg/(m3·s)或kJ/(m3·s)

T—— 温度，K

α—— 体积分数

β—— 时间松弛因子，s-1

ρ—— 密度，kg/m3

μ—— 黏度，Pa·s

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

下角标

cv—— 冷凝段流体

ev—— 蒸发段流体

E—— 能量

k—— 第k相

l—— 液相

M—— 质量

sat—— 饱和状态

v—— 气相

wc—— 冷凝段壁面

we—— 蒸发段壁面