高温高压氢气在三角凹槽及管内流动换热研究

刘 林，房玉良，武俊梅，王成龙，马 元，田文喜

(1.西安交通大学 航天航空学院 机械结构强度与振动国家重点实验室 陕西省先进飞行器服役 环境与控制重点实验室，陕西 西安 710049；2.西安交通大学 能源与动力工程学院 动力工程多相流国家重点实验室，陕西 西安 710049；3.西安航天动力研究所，陕西 西安 710100)

0 引言

随着人类空间探索规模的不断提高，对航天器的飞行时间、飞行载荷和飞行加速度等性能指标要求越来越高，因而对推进技术提出了新的要求。目前较为成熟的化学火箭推进系统受到燃料化学能和结构耐热性能限制，最大理论比冲约520 s，且能量密度也低，难以适应未来空间活动的动力需要。而电推进系统尽管其比冲可以高达上万秒，但它的推力水平却很低，也不能适应大规模深空探测任务[1-2]。相比之下，核热推进系统具有比冲高、推力大、控制灵活和工作时间长等特点，在深空探测任务中具有广阔的应用前景[3-4]。在核热推进系统中，工质在堆芯中吸收裂变热而被加热至高温高压状态，经尾喷管被加速到超音速流，向外喷出产生推力。氢气(临界压力1.296 4 MPa，临界温度33.145 K)因其分子量低，换热性能好，用作核热火箭的推进剂和堆芯冷却剂，不仅可以产生更大的比冲，同时可以简化系统。20世纪50年代始，美国和苏联着手研究核热推进技术，经过半个多世纪的不懈努力，固态堆芯核热推进技术得到一定程度的发展。固态堆芯设计中冷却通道有轴向、径向两种方案。NERVA型堆芯是轴向型的典型代表，堆芯燃料元件内分布着细长的轴向冷却剂通道[5]，作为冷却剂的氢气在细长的通道内流动并吸收堆芯的裂变热。

为获得更大的比冲，要尽可能提高堆芯的热流密度，并要求裂变热能够高效地传递给冷却剂[6]。堆芯热流密度大，冷却通道长，氢气压力高，流过堆芯后温升大，物性变化范围大，导致堆芯热工特性复杂，因此很有必要进行变物性条件下氢气在冷却通道内的流动换热特性研究。

氢气在高温高压(甚至超临界压力)条件下的实验难度高，公开文献报道也不多。McCarthy和Wolf[7]进行了氢气在圆管中的流动换热实验,实验进行的压力范围为0.22～9.32 MPa，壁面温度范围为461～1 245 K。Taylor开展了高热流密度条件下氢气在圆管中的流动换热实验，实验进行的压力范围为0.28～0.54 MPa，壁面温度范围为410～3 130 K[8]。Hess与Kunz进行了超临界氢气的流动换实验，压力范围为1.62～5.14 MPa，壁面温度范围为167～500 K[9]。基于以上实验结果，学者们得到了相应的实验关联式。关于氢气在真实堆芯通道中的流动换热实验公开报道很少，Lyon进行了氢气在Nuclear Furnace 1(NF—1)型核热推进系统堆芯冷却通道内的流动换热实验，冷却通道出口压力3.2 MPa，氢气出口温度高达2 450 K，得到了冷却通道壁面温度沿轴向的分布，但给出的实验数据很少[10]。

关于氢气在直通道内流动换热的数值模拟研究也有一些。Appel设计了一种核热推进堆芯燃料元件并对氢气在圆形冷却通道内流动换热特性进行了数值模拟研究，湍流模型采用Realizablek—ε模型，文中采用压力为3.5 MPa下的氢气物性数据，将氢气物性拟合成温度的单值函数[11]。Singh、Akyuzlu对圆通道内氢气的流动换热特性进行了研究，湍流模型采用标准k—ε模型，也考虑了氢气的物性随温度的变化，对影响圆通道内氢气流动换热的因素进行了分析[12-13]。艾青等研究了通道几何尺寸对氢气换热的影响[14]。房玉良等开展了高温、高流速氢气在圆管中的流动换热特性研究，采用SSTk—ω湍流模型，考虑了高温热解对氢气物性的影响，分析了进口流量和壁面热流密度对流动换热的影响[15]。Gould等研究了微细管道内高温氢气的流动换热特性，主要研究了三种管道(截面形状为正三角形、正方形和圆形)中温度场和速度场的分布情况[16]。Ji等对超临界氢气在长直圆管内的流动换热特性进行了数值研究，分析了传热恶化现象产生的原因[17-18]。

以上研究多针对圆管中氢气的流动换热特性，对如何改善氢气在冷却通道内综合流动换热性能的研究较少。本文开展了高温高压条件下氢气在带三角凹槽强化圆管内的流动换热数值模拟研究。

1 计算模型及方法简介

1.1 几何模型

本文首先对Lyon进行的NF—1型核热推进系统堆芯冷却通道内高温氢气换热实验进行了数值模拟，将计算结果与其实验结果进行对比，以检验数值计算模型和方法[10]。文献[10]中的冷却通道简化为直径D=2.3 mm，长度L=1 320 mm的长直圆管。进口质量流量为1.158 g/s，进口温度为340 K，出口压力为3.2 MPa。NF—1型核热推进系统堆芯冷却通道壁面热流密度沿轴向是变化的。在此基础上，本文模拟了沿轴向均匀热流加热的长直圆管和带三角凹槽强化圆管内氢气的流动换热情况，加热管段直径和长度与上面的相同，但为了使加热段达到速度充分发展，计算模型包括80 mm长的入口段。图1所示为长直圆管的几何结构示意图。带三角凹槽强化圆管分为内凹型和外凸型两种，沿轴向等距分布着槽宽w=1.15 mm的三角凹槽，如图2所示。主要分析了凹槽管相对于光滑圆管的强化换热能力以及凹槽管凹槽间距比I/D、凹槽深度比d/D、凹槽高度比H/D、凹槽顶点距离凹槽左边缘的间距比s/w对凹槽管强化换热能力的影响。I/D取4.35、8.70和28.70。在保持凹槽径向截面形状为等腰三角形的条件下，d/D和H/D取0.05、0.1和0.15。另外，保持d/D=0.05不变，s/w取0.25、0.5和0.75。

图1 长直圆管Fig.1 Long straight tube

图2 带三角凹槽强化圆管Fig.2 Triangular fluted enhanced tubes

1.2 控制方程及边界条件

本文采用轴对称坐标系下的二维稳态控制方程，忽略化学反应、辐射换热等，控制方程为

质量方程

(1)

动量方程

(2)

能量方程

(3)

式中：i、j=1、2，分别为轴向和径向的两个坐标方向；u为时均速度；u′为湍流脉动速度，m/s；T为温度；T′为湍流脉动温度，K；ρ为流体密度，kg/m3；μ为流体动力黏性系数，kg /(m·s)；cp为流体定压比热容， J/(kg·K)；λ为流体热导率，W/(m·K)。因为是湍流，动量方程中包括脉动引起的雷诺应力项，能量方程中包括脉动引起的附加热流，值得一提的是因为工质的高速流动和高温下黏性系数的增大，能量方程中还考虑了黏性耗散项Φ。

计算区域进口选用质量流量进口和均匀温度进口边界条件；出口选用定压力边界条件。所有算例的进口温度均为340 K，出口压力均为3.2 MPa；壁面为定热流密度且无滑移。

1.3 网格划分和网格无关性验证

长直圆管算例的长度为1 320 mm，直径为2.3 mm，验证工况为壁面热流密度q=2 MW/m2、进口质量流量m=1.158 g/s、进口温度Tin=340 K、出口压力pout=3.2 MPa。以出口主流总温Tb,total,out、进口静压pstatic,in以及流体的出口主流平均速度ub,out作为评判标准。计算结果如表1所示，在保证计算准确性同时考虑节约计算资源的情况下，使用35×7 000的网格模型。对于带三角凹槽强化圆管的局部网格划分如图3所示，通过网格加密处理，使得y+接近于1,以满足湍流模型需求。

表1 网格无关性验证Tab.1 Grid independence verification

图3 带三角凹槽强化圆管网格划分Fig.3 Meshing of triangular fluted enhanced tubes

本文采用Coupled算法进行压力与速度的耦合求解，采用二阶迎风格式进行空间离散。湍流模型采用SSTk-ω模型，收敛判定标准为残差小于10-7，其中能量项残差小于10-8，同时检测进出口流量守恒。

1.4 氢气物性模型

本文采用FLUENT内置的Real-Gas-Aungier-Redlich-Kwong 模型，具体如下

(4)

(5)

(6)

α(T)=α0(T/Tc)-n

(7)

n=0.498 6+1.173 5ω+0.475 4ω2

(8)

(9)

式中：p为压力，Pa；R为氢气的气体常数，4 124.3J/(kg·K)；v为比容，m3/kg；b、c、α0和n为系数；α(T)为温度修正因子；ω为氢气的偏心因子，-0.219；下标c为临界状态，Tc、pc和vc分别为33.145 K、1 296. 4 kPa和0.031 988 m3/kg。

氢气的热力学和输运性质受到温度和压力的共同影响。图4所示为基于NASA的氢气物性数据[19]所绘制的不同压力下的氢气物性随温度的变化曲线，可以看出压力对氢气的热导率、定压比热容的影响主要在高于2 000 K以后。对于本研究，高温氢气存在于冷却通道的出口区域，所以氢气的物性随温度变化的曲线拟合选取出口压力条件下的数据。如前所述，本计算的管道出口压力基于文献[10]的数值，大小为3.2 MPa，而3 MPa和3.5 MPa下氢气的物性相差很小，因此本文采用NASA的氢气物性数据[19]中压力为3.5 MPa、温度范围为300～3 000 K的氢气物性数据拟合了如下的氢气物性模型，在FLUENT的物性UDF模块中写入。

图4 氢气物性参数随温度的变化关系Fig.4 Temperature dependence of hydrogen properties

动力黏性系数

μ(T)=3.251 78×10-6+2.068 54×10-8T-

4.149 57×10-12T2+5.558 3×10-16T3

(10)

热导率

λ(T)=0.012 13+7.672 66×10-4T-

7.891 87×10-7T2+7.074 22×10-10T3-

3.064 47×10-13T4+5.215 54×10-17T5

(11)

定压比热容

cp(T)=14 618.54-0.482 12T-5.513 41×

10-4T2+2.832 85×10-6T3-1.757 74×

10-9T4+3.520 04×10-13T5

(12)

1.5 数据处理

本文各特征参数定义如下：

雷诺数

(13)

表面传热系数(局部)

(14)

努赛尔数(局部)

(15)

摩擦因子

(16)

Δp=pin-pout=Δpa+Δpf

(17)

(18)

综合强化传热因子[20]

(19)

式中：D为冷却通道的直径，m；q为加热面的热流密度，W/m2；h为表面传热系数，W/(m2·K)；Tw为壁面温度，K；Tb为流体平均温度，K；Δp、Δpf和Δpa为进出口压降、摩擦压降和加速压降，Pa；L为加热段的总长，m；m为流体的质量流量，kg/s；A为冷却通道截面面积，m2；下标in、out、ave和0分别代表进口、出口、进出口平均值和光滑圆管。由于冷却通道内氢气的密度变化剧烈，本文计算摩擦因子时考虑了加速压降Δpa。

2 结果与讨论

2.1 模型与方法验证

文献[10]中模拟了NF—1型核热推进系统堆芯裂变热释放特征，沿冷却通道轴向的热流密度是变化的，如图5所示。实验结果给出了通道壁面温度沿轴向的变化，图6展示了本文采用两种不同的湍流模型所计算出的壁面温度与文献[10]实验结果的对比情况，可以发现，对于光滑通道而言，两种湍流模型的计算结果很接近。本文的研究内容是带三角凹槽强化圆管内的对流换热特性，SSTk—ω模型对于壁面凹槽附近的细微流场特征的捕捉要优于k—ε模型，因此后面的计算选用SSTk—ω湍流模型。从图6可以看出，计算结果和该实验的最大相对误差为15%，壁面温度计算值高于实验值，偏差随着壁面温度的上升而增大。出现偏差的主要原因是计算中没有考虑壁面包壳材料的热阻，忽略了固体壁面的导热和辐射换热，这两项的影响随着壁面温度升高而增大。总体来说，该计算结果是可以接受的，采用本数值计算方法可以进行高温高压条件下氢气的流动换热特性研究。

图5 热流密度分布[10] Fig.5 Heat flux distribution[10]

图6 不同湍流模型的计算值和实验的壁温对比Fig.6 Comparison of the wall temperature distribution between calculated values and experimental data for different turbulence models

2.2 冷却通道内氢气的流动换热特性研究

图7为计算得到的氢气在NF—1型核热推进系统堆芯冷却通道内的温度场和马赫数分布云图，可以看出，氢气的温度沿轴向不断上升，且沿轴向壁面温度的升高速率大于通道中心区域氢气温度的升高速率。这是因为氢气在高速、高温条件下流动时，动力黏性系数很大，因粘性产生的气动加热效应很强，使得壁面附近气体温度和壁面温度快速升高。计算结果显示在文献[10]的冷却通道内，氢气的马赫数沿轴向不断上升，这是因为氢气被加热，温度不断升高，氢气发生膨胀，氢气的流速不断提高，但同时因氢气温度不断升高，压力也相对较高，当地声速沿轴向也不断提高，所以马赫数不断增大，但始终处于亚声速流动，近壁面区由于黏性作用流速较低。

图7 氢气温度场和流场分布云图Fig.7 Contours of hydrogen temperature and flow field distribution

图8展示了局部努赛尔数沿轴向的分布情况，冷却通道进口处由于边界层厚度较小，换热较强，局部努赛尔数较大，之后因为边界层发展过程中厚度逐渐增大，且随着氢气和壁面不断地进行热量交换，氢气温度不断升高，进而氢气的热导率不断增大，局部努赛尔数不断减小。

图8 局部努赛尔数轴向分布Fig.8 Axial distribution of local Nusselt number

2.3 带三角凹槽强化圆管内氢气的流动换热分析

对于带三角凹槽强化圆管内氢气流动换热的数值模拟，保持壁面热流密度q=2 MW/m2不变，在进口温度Tin=340 K，出口压力pout=3.2 MPa条件下，将进口质量流量m=1.158 g/s作为基础工况，质量流量按±10%、±20%进行变化。同时计算了相同直径的光滑圆管在不同质量流量(Re数)下的平均努赛尔数Nu0用于比较。Nu/Nu0表示带三角凹槽强化圆管和相同直径的光滑圆管的平均努赛尔数的比值，f/f0表示带三角凹槽强化圆管和相同直径的光滑圆管的摩擦因子的比值。为表示带三角凹槽强化圆管相对于相同直径的光滑圆管的综合流动换热特性，引入了综合强化传热因子PEC。

2.3.1 凹槽深度/高度对流动换热的影响

图9所示为进口质量流量m=1.158 g/s，凹槽间距比I/D=8.70，凹槽宽度比w/D=0.5，凹槽深度比d/D和凹槽高度比H/D取0.05、0.1和0.15时凹槽附近的局部流线图和温度云图。可以看出，氢气流经内凹型三角凹槽时，在凹槽顶点后形成横向涡，凹槽深度越大，横向涡越强，横向涡引起顶点后壁面当地换热增强，所以壁面温度降低。而且内凹型带三角凹槽强化圆管内沿轴向持续的周期性凸起对主流产生周期性扰动，使主流流线弯曲。氢气流经外凸型三角凹槽时，槽内形成横向涡，凹槽高度越大，横向涡越强，但因为横向涡几乎封闭在槽内，对主流扰动很小，主流流线沿轴向几乎平行。该处通过壁面进入的热量需要通过凹槽内的氢气再传递给主流区的氢气，使得热边界层变厚，当地的加热量不能被及时带走，壁面及壁面附近流体温度升高。这说明外凸型带三角凹槽强化圆管对主流的扰动较小，强化换热能力有限。

图9 不同凹槽深度/高度的带三角凹槽强化 圆管的局部流线图和温度云图Fig.9 Local contours of temperature and streamlines of triangular fluted enhanced tubes with different groove depth or height

图10所示为3种不同凹槽深度的内凹型带三角凹槽强化圆管的传热特性曲线。结果表明，内凹型带三角凹槽强化圆管的强化换热能力随进口质量流量(Re数)的增大变化很小。但随着槽深的增大，内凹型带三角凹槽强化圆管的强化换热能力增强。原因在于如图9所呈现的在三角凹槽顶点后面的横向涡的产生，周期性的涡使得流体扰动加剧，同时破坏了边界层的发展，强化了换热。相对于光滑圆管，在凹槽深度比d/D=0.15时，内凹型带三角凹槽强化圆管的努赛尔数提高达53%，且凹槽深度越大，上述影响越剧烈，换热效果越好。

图10 不同凹槽深度的内凹型带三角凹槽 强化圆管的传热特性变化曲线Fig.10 Nu/Nu0 of inward triangular fluted enhanced tubes with different groove depth

图11所示为3种不同凹槽高度的外凸型带三角凹槽强化圆管的强化换热能力变化曲线。结果表明，外凸型带三角凹槽强化圆管的强化换热能力随进口质量流量(Re数)的增大变化也很小。随着槽高的增大，外凸型带三角凹槽强化圆管的努赛尔数先增大后减小。凹槽高度比H/D=0.1时，外凸型带三角凹槽强化圆管的强化换热能力最强，但Nu/Nu0也小于1.041，强化换热效果不显著。这也是与图9所呈现出的流场相关，高度小的凹槽对主流产生些许扰动，对强化换热有利，但努赛尔数的提升幅度很小。随着凹槽高度的增加，流体在凹槽内的扰动会增强，但随着槽高的进一步增大，滞留在凹槽内的氢气主要通过热传导的方式与主流氢气进行热量交换，壁面对流换热能力减弱。

图11 不同凹槽高度的外凸型带三角凹槽强化圆管 的传热特性变化曲线Fig.11 Nu/Nu0 of outward triangular fluted enhanced tubes with different groove height

图12为3种不同凹槽深度的内凹型带三角凹槽强化圆管的阻力特性变化曲线，结果表明，随着进口质量流量(Re数)的提高，相对于光滑圆管，三种凹槽管的阻力都不同程度地提升，且槽深越大，提升越明显。阻力增大的原因在于内凹型带三角凹槽强化圆管内横向涡的产生，而横向涡的强度与槽深和氢气主流流速相关，凹槽越深、流速越大，横向涡越强，阻力提高越显著。凹槽深度比d/D=0.15，在Re=80 000时，内凹型带三角凹槽强化圆管的摩擦因子约为光滑圆管的14倍。

图12 不同凹槽深度的内凹型带三角凹槽 强化圆管的阻力特性变化曲线Fig.12 f/f0 of inward triangular fluted enhanced tubes with different groove depth

图13为3种不同凹槽高度的外凸型带三角凹槽强化圆管的阻力特性变化曲线，其变化趋势和内凹型一致，但摩擦因子增加幅度较小。凹槽高度比H/D=0.15，阻力提升约11%。同样是随着凹槽高度的增大，其阻力增加幅度也变大。前面的研究结果表明，外凸型带三角凹槽强化圆管的强化换热能力有限，本文暂不深入研究。而内凹型带三角凹槽强化圆管可以带来显著的强化换热效果，但也带来了阻力系数的大幅增加。为了掌握内凹型带三角凹槽强化圆管的综合流动换热特性，对不同管道在不同进口流量下的综合强化传热因子PEC进行了比较。图14所示为3种不同凹槽深度的内凹型带三角凹槽强化圆管的综合换热性能变化曲线，结果表明，随着进口雷诺数的增大，其综合强化传热因子变小，且槽深越大，综合强化传热因子越小。

图13 不同凹槽高度的外凸型带三角凹槽强化 圆管的阻力特性变化曲线Fig.13 f/f0 of outward triangular fluted enhanced tubes with different groove height

图14 不同凹槽深度的内凹型带三角凹槽强化 圆管的综合换热性能变化曲线Fig.14 Overall heat transfer performance comparison of inward triangular fluted enhanced tubes with different groove depth

2.3.2 凹槽间距对流动换热的影响

保持凹槽宽度比w/D=0.5，凹槽深度比d/D=0.05不变，以内凹型带三角凹槽强化圆管为例，研究凹槽间距对流动换热的影响。图15所示为三种不同凹槽间距的内凹型带三角凹槽强化圆管的综合换热性能变化曲线。结果表明，随着进口雷诺数的增大，其综合强化传热因子不断减小，凹槽间距越小，综合换热性能越弱。原因在于凹槽间距越小，凹槽数目就越多，换热增强的同时阻力明显增大。

图15 不同凹槽间距的内凹型带三角凹槽强化圆管 的综合换热性能变化曲线Fig.15 Overall heat transfer performance comparison of inward triangular fluted enhanced tubes with different groove spacing

2.3.3 凹槽顶点位置对流动换热的影响

保持凹槽宽度比w/D=0.5，凹槽深度比d/D=0.05和凹槽间距比I/D=8.70不变，以内凹型带三角凹槽强化圆管为例，研究凹槽顶点距离凹槽左边缘的间距比s/w对流动换热的影响。图16所示为进口质量流量m=1.158 g/s时，不同凹槽顶点位置的内凹型带三角凹槽强化圆管的局部流线图和温度云图。结果表明，相较于凹槽顶点位于中点时(图16(b))，凹槽顶点左移(图16(a))时氢气接触凹槽顶点时流速最低，在凹槽后形成的横向涡能够很好地覆盖凹槽后壁面，对换热有利，另外氢气流经凹槽时与主流方向产生的夹角最大，使得阻力提升。凹槽顶点右移(图16(c))时横向涡变强，但由于氢气流经凹槽时流速提升，横向涡后移，使得凹槽后局部壁温升高，流速提升也会导致凹槽对氢气的阻力增大。

图16 不同凹槽顶点位置的内凹型带三角凹槽强化 圆管的局部流线图和温度云图Fig.16 Local contours of temperature and streamlines of inward triangular fluted enhanced tubes with different vertex position

图17所示为3种不同凹槽顶点位置的内凹型带三角凹槽强化圆管的综合换热性能变化曲线，不难发现，相较于凹槽径向截面形状为等腰三角形也就是凹槽顶点位于中点时，凹槽顶点的左右移动使得综合强化传热因子减小，且凹槽顶点右移时综合换热性能最差。

图17 不同凹槽顶点位置的内凹型带三角凹槽 强化圆管的综合换热性能变化曲线Fig.17 Overall heat transfer performance comparison of inward triangular fluted enhanced tubes with different vertex position

3 结论

1)本文采用Coupled算法、SSTk-ω湍流模型，将氢气的物性看作温度的单值函数进行高温高压条件下氢气在管内的流动换热数值模拟方法可行。

2)NF-1型核热推进系统堆芯裂变热沿冷却通道轴向大致呈抛物线变化，通道内的氢气流动属于亚音速湍流，冷却通道壁面温度沿轴向先逐渐上升在靠近出口时略微下降。

3)内凹型带三角凹槽强化圆管随槽深的增大换热增强，阻力显著提升。外凸型带三角凹槽强化圆管随槽高的增大换热先增强后减弱，阻力提升不大。内凹型带三角凹槽强化圆管的强化换热能力优于外凸型带三角凹槽强化圆管，且随槽深的增大，综合换热性能降低。

4)凹槽间距的减小有利于强化换热，但综合换热性能降低。非对称三角凹槽的综合换热性能低于对称三角凹槽，其中凹槽顶点右移时综合换热性能最差。