高温两相闭式热虹吸管的瞬态气液相变数值模拟研究

2022-07-11耿伟轩姚寿广

江苏科技大学学报(自然科学版) 2022年2期

沈妍，刘锐*，耿伟轩，于萍，姚寿广

(1.江苏科技大学机电与动力工程学院张家港 215600) (2.江苏省秦淮河水利工程管理处南京 210000) (3. 徐州工程学院机电工程学院徐州 221018) (4.江苏科技大学能源与动力学院镇江 212100)

高温热管常采用钠、钾、钠钾混合等碱金属工质[1]，工作温度为450 ℃以上，具有高效等温传热能力[2].高温热管由空间技术应用衍生到石油化工、新能源利用、冶金工业等领域，已有50多年的研究发展历史.高温两相闭式热虹吸管内无吸液芯，靠重力回流促进气液循环[3].针对两相闭式热虹吸管的各种传热性能实验研究[4-8]较多，在运行时，涉及封闭空间内钠工质的融化、蒸发、冷凝和流动等过程[9]，且各过程互相影响，其内部传热传质机理模拟是一个研究难点.

气液两相流可视化数值模拟的发展，使热管内气液相变传热模拟成为可能.采气液可视化相变数值模拟模型主要有流体体积法(volume of fluid，VOF)模型[10]、混合模型[11]和欧拉模型.其中，VOF模型较混合模型能很好模拟热管的工作状态[12]，得到清晰的气液两相界面.相变传质模型主要有Sharp interface 模型，Hertz Knudsen Schrage相变模型[13]，Lee相变模型[14]和不凝性气体壁面冷凝模型等.研究者采用VOF模型和相变传热模型对两相闭式热虹吸管进行了模拟研究，文献[15]考虑了水重力热管中冷凝段存在不凝性气体对传热性能的影响；文献[16]分析不同加热功率下重力热管的传热性能，结果显示在相同加热功率下，实验数据和模拟结果具有高度一致性；文献[17]对脉冲加热功率下的热虹吸管热响应进行了研究，得到了充液率为30%时热管冷凝段的液膜分布情况和不同时刻轴向温度分布，分析了不同脉冲热流密度变化时的热管稳定时间；文献[18]对制冷剂R134a 和 R404a为工质的重力热管热力性能进行模拟，得到热管内不同制冷剂或纯水作为换热工质不同的沸腾现象；文献[19]研究了Lee模型中蒸发系数分别为1、0.95、0.9和0.8时的热管传热性能，模拟结果与实验值对比，发现蒸发系数为0.9时吻合最好.此外，通过模拟结果的壁温分析得到随着功率100～200 W增大时，蒸发段和冷凝段热阻下降，且比较了有无考虑瞬时时间松弛参数的模拟结果，考虑了瞬时时间松弛参数的模拟结果与实验值误差更小.文献[20]对热管内插入翅片的热虹吸管进行了数值模拟，分别采用VOF模型和集总参数模型进行研究，并进行了场协同的强化传热分析；文献[21]对不同倾角铜-水重力热管性能影响数值研究；文献[22-23]对两相封闭式水热虹吸管中间歇沸腾传热进行模拟，引入了初始过热度，又对用于浅层地热利用的氨热管进行了可视化数值模拟和实验研究[24].现有研究基于CFD软件Fluent，采用VOF模型和UDF编写蒸发冷凝模型，主要进行了水、R134a、氨、纳米流体[25]等工质的中低温两相闭式热虹吸管的数值模拟.

高温热管采用碱金属钠为工质，由于液钠具有高导热系数，高沸点的优良物性，其初始过热度、临界热流密度和沸腾滞后等是钠沸腾机理研究的重点.钠有强腐蚀性，与氧气会发生化学反应，对高温热管制造密封性和运行稳定性提出更高的要求.热管内钠的气液相变、温度分布对高温热管的失效有重要影响，液钠在密闭的高温热管内沸腾传热现象模拟还鲜有报道，液钠的物性模型，和质蒸发冷凝模型还不清楚.文中建立了高温热虹吸管的二维模型进行可视化数值模拟研究，分析高温热管内气液相变现象和规律，为高温热管的模拟研究奠定基础.

1 模型建立

1.1 物理模型

高温热虹吸管总长为1 m，外径为0.019 m，壁厚0.75 mm.蒸发段长为0.4 m，绝热段长为0.2 m和冷凝段长为0.4 m.

1.2 数学模型

数学模型假设：

(1) 流动域中层流流动.

(2) 工质的蒸汽态为不可压缩理想气体，采用理想气体状态方程计算.

采用连续性方程、动量方程，能量方程、层流模型、流体体积分数模型(VOF)、连续表面力模型对重力热管进行数值计算.

连续性方程为：

(1)

式中：ρ为密度；t为时间；u为速度矢量；为向量微分算子；SM为有蒸发或冷凝引起的质量源相.

增加含气相体积分数αv的连续性方程为：

(2)

模拟过程中只涉及气液相变，两相体积分数和为1.

αl+αv=1

(3)

式中：αv和αl分别为气相体积分数和液相体积分数.

主相液相的体积分数αl根据式(3)求解.

动量方程为：

ρg+FCSF

(4)

式中：μ为流体黏度；ρ为流体受到的压强，Pa；g为重力加速度；FCSF是表面张力体积力，Pa/m.动量方程中添加源项来考虑表面张力的影响.表面张力用表面的压力跳跃来表示.使用散度定理，将表面上的力作为体积力获得.该体积力FCSF采用连续表面力模型[26]模拟为：

(5)

能量方程为：

(6)

el=Cp,l(T-Ts)

(7)

ev=Cp,v(T-Ts)

(8)

式中：e为单位质量流体的能量；T为流体温度；k为有效导热系数；SE为能量源相；Cp为比热；Ts为饱和温度.

动量方程和能量方程中混合相物性为：

ρ=αvρv+(1-αv)ρl

(9)

μ=αvμv+(1-αv)μl

(10)

k=αvkv+(1-αv)kl

(11)

(12)

气液相变模型采用在连续性方程和能量方程中添加源项实现，质量方程源项分为液相源项和气相源项.采用与Lee传热传质模型相似的形式，假设相变在压力不变的准静态平衡饱和温度下进行.

液相源项：

蒸发时，

(13)

冷凝时，

(14)

气相源项：

蒸发时，

(15)

冷凝时，

(16)

液相和气相能量方程源项分别为：

(17)

(18)

式中：β1为蒸发系数；β2为冷凝系数；Tmix为混合流体温度，Lfg为气化潜热.

1.3 工质物性

采用1 000 K时钠工质液态物性参数，密度为778 kg/m3，定压比热容Cp为1 260 J/(kg·K)，导热系数56.6 W/(m·K)，粘度1.809×10-4Pa·s.液态金属钠在沸点1 000 K时汽化潜热为4 026 kJ/kg.

1.4 边界条件和模型设置

当液钠在316L不锈钢260 ℃时，接触角为140°，为非润湿工况.当温度超过260 ℃时，接触角迅速下降.当300 ℃时，液钠在不绣钢表面的接触角为90°，为润湿工况[27].壁面上的液钠接触角随着体系温度、传热传质过程、平衡时间等因素变化着.考虑管内壁面粗糙度的影响，文中取316L表面470 ℃时接触角60°作为热管内液钠接触角设置值.沸点为1 000 K，充液率为0.5.蒸发段外壁面热流密度为83 766 W/m2.当蒸发系数为0.1，冷凝系数为0.2，也可根据网格尺度和时间步长，通过尝试性计算后确定系数[28].热管从常温到高温加热过程模拟缓慢，且内部气液相变过程主要在达到相变温度后发生.为了减少模拟时间，不考虑升温过程，初始温度设为999 K，得到简化启动加热升温过程后的高温热管蒸发段瞬时沸腾现象.液体充液率为0.5，即轴向0～0.2 m内为液体工质，0～0.4 m内为受热的蒸发段，0.4～0.6 m内为绝热段，0.6～1 m之间为对流换热的冷凝段.冷凝段对流传热系数为1 478 W/(m2·K).环境温度为298 K，初始压力为0，库朗数Courant数设置为0.5.

计算时，首先在气体体积分数αv>0的单元，计算αv的梯度，根据梯度判断气体分布属于20种类型的哪一种；根据相界面类型求解式(2)，通过计算单元界面的流量，并更新下一个时刻每个单元的气体体积分数αv；根据气液流体体积分数，通过式(9～12)对密度和粘度等物性参数进行加权计算；计算与流体体积分数相关的受力，连续表面力式(5)，获得连续性方程、动量方程和能量方程与单相流动情况一致，按单相流动的求解方案对方程进行求解.文中采用几何重建法对位于两相之间界面附近的单元格进行界面处理.压力速度耦合采用SIMPLE方法.压力离散采用PRESTO！方法，动量和能量离散采用一阶迎风格式[15-16].时间步长0.000 1 s.

蒸汽的雷诺数Rev为：

(19)

式中：rhv为蒸汽腔水力半径，取0.008 75 m；Q为输入热量，取1 999 W；Av为蒸汽流道截面积，取2.4×10-4m2；μv为蒸汽黏度，即2.211×10-5Pa·s；Lfg为汽化潜热4 026 kJ/kg.计算得蒸汽的Rev为1 637>2 300，为层流流动.

2 数值模拟结果

2.1 网格无关性分析和模型验证

根据不同单元数的外壁面上各段的平均温度，计算了不同网格数热阻如图2，热阻变化在0.176～0.193之间，随着网格数增加，热阻趋于稳定.由于Courant数与时间步长和网格尺寸有关，若网格越细，Courant数更大，收敛速度逐渐加快，但稳定性更差，越难收敛，所需时间步长更小，计算时间增长.

图2 不同网格数的热阻

为了使Courant数保持较小范围，没有选用满足网格无关性条件的最细网格数，选用了次之的78 430网格数，减少计算时间成本.

模拟值与实验值[29]比较如图3，最大偏差为17.8%，主要是蒸发段偏差较大，实验中的管壁温度散热损失等较大，蒸发段温度偏低.

图3 热管模型模拟计算值与实验值比较

2.2 高温热虹吸管瞬时气液相变模拟分析

2.2.1 钠热虹吸管沸腾

高温热管沸腾气化现象如图4.

图4 高温热管蒸发段瞬时变化

蒸发段壁面处先形成气化核心，产生蒸汽后扰动液池，气泡合并溢出、壁面处原气泡位置由液体补入.经换算，液钠工质升高1 K所需的加热时间为0.051 2 s，因此，在0.05 s时还未见液体蒸发气化现象.气化所需一定过热度，0.1 s时仍未见气化.加热后0.135 s时靠近壁面的液体受热先达到气化温度，此时为饱和沸腾，说明达到一定过热度后发生相变现象.

热管加热时不同时刻轴心速度分布如图5.当0.05 s时，热管内轴中心线上流速几乎为0，主要使液体加热升温，还没有产生蒸汽.当加热0.2、0.4 s时，在0～0.2 m的工质充液段内接近0且变化不大，充液段壁面附近处工质初步受热产生气泡，轴心液体流速较低.在0.2～0.6 m的蒸汽段，热管轴心气流速度随着加热时间逐渐增大，且波动剧烈.由于管内充斥着低压蒸汽，产生的少量蒸气受热上升，流速不高，易受到壁面的摩擦形成环流，因此在轴心速度上呈较大波动.在0.6～1 m的冷凝段壁面对流换热使蒸汽降温，速度波动且下降.在靠近冷凝段顶端1m处，气流开始冷凝产生局部低压使蒸气流速剧烈增加再冲击管顶，流速降为0.当加热0.6 s时，工质轴心处速度相对稳定，尤其是蒸汽段速度无波动现象.这是由于随着加热的进行，蒸汽量逐渐增大，形成稳定的气流.

图5 不同时刻轴心速度分布

工质在0～0.15 m充液段内液体的轴心速度逐渐变大，在0.15～0.3 m蒸发段与绝热段交界处，气液混合物流速波动剧烈并逐渐上升，在接近蒸发段尾端0.37 m处气流速度达到最大，出现了蒸发段对工质的加速现象，在绝热段气流速度基本不变，在0.6～1.0 m冷凝段内轴心气流速度逐渐下降，体现了冷凝段内工质冷凝，蒸汽量减少对气流减速的作用.

2.2.2 钠热虹吸管冷凝

当热管吸热蒸发后，冷凝液体在壁面形成，经计算冲液量为蒸发段体积的0.5时，若液膜遍布热管内壁面时，液膜平均厚度为1.75 mm，由于热管蒸发段有时有少量液体，冷凝液膜有厚有薄，有的内壁面无液膜，有的形成液滴，如图6.

图6 热管稳定运行的气体体积分数分布

(1) 冷凝段

冷凝段气液相变分布和流线分布如图7.为了显示清晰，图中未显示壁面固体壁面部分，只显示流体部分.图7(a)可以看出，气流从蒸汽腔中间上升到冷凝段顶端向两边散开在壁面上冷凝形成液膜回流.图7(b)可以看出，在冷凝段顶端，液膜逐渐形成，随着回流的进行，壁面上冷凝液膜逐渐流下，液膜先变厚再变薄，偶尔有液膜断层和液膜厚度变化，最厚的地方液膜约1 mm.图7(c)可以看出，在靠近绝热段的冷凝段，壁面附近处明显有气液逆流现象.当气流速度较大时，发生液体夹带，使冷凝液离开壁面形成含气率0.35～0.8的小液滴.根据蒸汽腔内流线，气流并不是直线上升，受到局部小液滴和液膜相界面表面力、壁面处液钠向下流动等因素的影响，有涡流形成，气体曲折上升，夹带液滴.

图7 冷凝段的气体体积分数和流线分布

(2) 绝热段

绝热段气液分布如图8，图8(a)可以看出，靠近冷凝段的绝热段0.55～0.6 m处，有液体被夹带向上的现象，液滴较多.图8(b)可以看出，靠近蒸发段的绝热段0.4～0.45 m处，无液滴，蒸汽从蒸发段出来加速到绝热段，气流上升速度大.

图8 绝热段气体体积分数和流线分布

(3) 蒸发段

蒸发段气液分布情况如图9.图9(a)可以看出，仅在0.395 m高度右侧壁面和0.37 m高度左侧壁面上有长度1 cm液膜，液膜很少，其余部分多为小液滴.图9(b)可以看出，液体回流后多在壁面上直接蒸发成气体.蒸发段上升气流流线均匀，涡流较少.图9(c)可以看出，蒸发段底部有少量液体留存.

图9 蒸发段气体体积分数和流线分布

2.3 高温热虹吸管瞬时传热性能分析

当热管初始温度为999 K时，加热后轴向上蒸汽温度和管子内壁面温度如图10.

图10 热管轴向温度随时间的变化

由于冷凝段与外界对流传热的作用，冷凝段内壁温度和轴心温度均下降，且内壁温度低于同时刻轴心蒸汽温度，温差约20 K.当0.629 s时，轴心蒸汽温度达到1 000 K，与蒸发段温度持平，冷凝段内壁面温度为840 K，温差为160 K，此时认为热管进入正常运行工况.

不同位置的液相体积分数如图11，当x=0.000 75 m位置位于热管左侧壁面，液相分布最多，在绝热段和冷凝段都有间歇的液膜形成，在蒸发段底端有少量液体存留，壁面处有不饱和液体存在.在x=0.001 75 m和x=0.002 75 m位置，即分别离壁面1 mm和2 mm处，液体分布明显减少，说明仅有部分位置的液膜达到了1 mm厚，且离壁面的位置越远，液体冷凝越少.在x=0.009 5 m位置即轴心处液相分布较少，轴心处的液体主要是个别漂浮的液滴.

图11 不同位置的液相体积分数分布

轴心处速度分布如图12.

图12 蒸汽腔轴线上速度分布

蒸发段和绝热段气流速度较大，在冷凝段内气流速度沿轴向减小.由于蒸发段气化现象扰动，轴心蒸汽速度不稳定，液滴蒸发使气流速度较大，而残存的液滴使该位置上速度较小.沿轴向蒸发段气流速度增大，绝热段气流速度不变，冷凝段速度逐渐下降.在热管运行中绝热段0.485 m位置处气流达到最大速度110.53 m/s，未达到声速.