本文首先对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

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

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

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

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

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

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

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

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

式中:p为压力,Pa; R为氢气的气体常数,4 124.3J/(kg·K); v为比容,m³/kg; b、c、α₀和n为系数; α(T)为温度修正因子; ω为氢气的偏心因子,-0.219; 下标c为临界状态,T_c、p_c和v_c分别为33.145 K、1 296. 4 kPa和0.031 988 m³/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^-12T²+5.558 3×10^-16T³(10)

热导率

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

7.891 87×10^-7T²+7.074 22×10^-10T³-

3.064 47×10^-13T⁴+5.215 54×10^-17T⁵(11)

定压比热容

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

10^-4T²+2.832 85×10^-6T³-1.757 74×

10^-9T⁴+3.520 04×10^-13T⁵(12)

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

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

文献[10]中模拟了NF—1型核热推进系统堆芯裂变热释放特征,沿冷却通道轴向的热流密度是变化的,如图5所示。实验结果给出了通道壁面温度沿轴向的变化,图6展示了本文采用两种不同的湍流模型所计算出的壁面温度与文献[10]实验结果的对比情况,可以发现,对于光滑通道而言,两种湍流模型的计算结果很接近。本文的研究内容是带三角凹槽强化圆管内的对流换热特性,SST k—ω模型对于壁面凹槽附近的细微流场特征的捕捉要优于k—ε模型,因此后面的计算选用SST k—ω湍流模型。从图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

图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

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