1. 引言

深空探索对于拓展人类生存空间具有重要意义。与化学电池和太阳能光伏阵列相比，空间核反应堆电源具有能量密度高、工作时间长、功率覆盖范围广、不受太空环境影响等特点，是一种很有前途的空间探索技术。目前，利用氦氙(He-Xe)气体混合物作为冷却剂的小型气冷反应堆继续引起全世界的关注。从普罗米修斯计划开始，美国建立了一个全面的He-Xe气冷空间反应堆计划。俄罗斯在2009年提出TEM项目，目标是开发兆瓦级He-Xe气冷反应堆。国内Qin [1]、Zhang [2]、Meng [3]、Ning [4]等人也提出了小型He-Xe反应堆的概念设计方案。

He-Xe反应堆中常用40 g/mol的He-Xe为冷却剂，其普朗特为0.21，比传统工质小(空气0.7)，同时比液态金属大很多(≈0.01)。普朗特数决定了流动和热边界层相对大小，会影响湍流换热特性。当普朗特数降低时，流动边界层厚度会逐渐小于热边界层，经典的雷诺比拟方法会发生失效，有研究表明，传统的传热关系式对低普朗特气体的对流换热预测偏差较大[5]。为了研究低普朗特数氦氙混合气体的流动传热特性，Taylor等人[5]对空气以及二元混合气体(He-Xe, H₂-Xe)在圆管内的传热特性进行了实验研究，提出了适用于计算氦氙混合气体的局部换热关系式。Nakoryakov [6]和Vitovsky [7]等人对40 g/mol氦氙混合气体在等边三角形通道中的传热特性进行了研究，得到的实验数据与Kirov等人[8]提出的对流换热关系式符合较好。Elistratov等人[9]实验研究了长圆柱形管中的传热，得到了传热系数随雷诺数的变化趋势。Qin等人[10]的圆管实验研究了He-Xe混合气体的对流换热特性，得到了新的关系式来预测努塞尔数。Vitovsky [11]实验研究了氦氙混合气体在加热的薄壁准三角形管内流动时的换热过程，提出了一种计算非定常流场局部换热系数的方法。He-Xe气体流动换热特性实验研究主要是基于简单通道内，但其研究成果为数值模拟研究提供了宝贵的实验数据基础。

在新型高效动力反应堆中，常采用紧密的燃料棒束布置设计，这种紧密栅的结构能提高堆芯功率密度，减小堆芯体积。因此，氦氙空间反应堆中紧密栅通道成为了不错的选择，而国内学者在设计堆芯结构时也有用到紧密栅结构[3]。因此，研究氦氙在紧密栅通道中的流动换热是有必要的。对氦氙混合气体的数值研究中，对圆管[12]和棒束子通道[13]的研究较多，但目前缺少紧密栅棒束通道结构对氦氙混合气体流动换热影响的研究，如不同栅径比结构及其强化换热的研究。对氦氙混合气体而言，现有的实验研究缺少对多棒束通道的实验，从而使实验成本昂贵。选择合适的湍流模型以及物性方法对He-Xe紧密栅结构中的流动传热进行CFD数值模拟能让我们更好地了解He-Xe在堆芯结构中的流动传热特性。

本文针对19棒紧密栅棒束通道中He-Xe混合气体的流动换热特性进行研究，通过建模与CFD数值求解进行分析研究，分析紧密栅结构中的流动换热情况，分析不同栅径比对He-Xe混合气体流动换热特性的影响。

2. 计算方法

2.1. CFD模型

在具有粘性的流动传热计算中，常用的湍流模型有k-ε模型、k-ω模型与SST k-ω模型。在各种湍流模型中，有研究表明，SST k-ω模型的He-Xe计算结果与实验符合较好，有较好的准确性[14] [15]。在本文的后续部分，证明了雷诺时均SST k-ω模型的准确性，并将其应用于本研究。SST k-ω控制方程如下：

$\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$ (1)

$\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_i}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\mu +{\mu }_t\right)\frac{\partial u_i}{\partial x_j}\right]+\rho g_i$ (2)

$\frac{\partial \left(\rho u_{j}h\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\lambda +\frac{{\mu }_t}{P{r}_t}\right)\frac{\partial T}{\partial x_j}\right]+{\tau }_{ij}{u}_i$ (3)

$\frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k-Y_k+S_k+G_b$ (4)

$\frac{\partial \left(\rho \omega u_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial \omega }{\partial x_j}\right]+G_{\omega }-Y_{\omega }+D_{\omega }+S_{\omega }+G_{\omega b}$ (5)

此外，由于He-Xe的低普朗特数特性，对SST k-ω模型中的湍流普朗特数进行了修改，采用Zhou等人[12]提出的Prt模型，如下：

$Pr={\left\{\frac{1}{2P{r}_{t\infty }}+\frac{a}{P{r}_{t\infty }^{0.5}}-a^2\left[1-\exp \left(\frac{-1}{aP{r}_{t\infty }^{0.5}}\right)\right]\right\}}^{-1}$ (6)

$a=\frac{0.3Pr{\mu }_t}{\mu }$ (7)

$P{r}_{t\infty }=0.8+\frac{30}{PrR{e}^{0.888}}$ (8)

2.2. 计算模型与网格划分

所建立的紧密栅带绕丝棒束通道如图1所示。通道内燃料棒的直径为10 mm，包壳厚度为1 mm，长度为600 mm；燃料棒呈三角形排列，间距为12 mm。网格的生成分成两个区域：固体区域和流体区域，对流体域壁面进行加密，确保大多数近壁面单元y⁺ < 1，网格质量在0.7以上，符合计算的要求。为确保网格数量对计算结果无影响，采用4组网格进行网格独立性验证，网格数量分别为414万、597万、804万、1050万。图2结果表明，网格在达到414万后，计算结果与网格数量无关，最后选用414万网格的模型进行后续计算。

Figure 1. Geometric model and mesh division of compact rod bundle channel

图1. 紧密栅棒束通道几何模型以及网格划分情况

Figure 2. Grid independence verification

图2. 网格无关性验证

2.3. 计算方法验证

为了验证湍流模型和物性计算方法的准确性，对2组圆管实验进行了对比验证。对比的实验为Taylor圆管实验[5]，计算模型如图3所示。在Taylor圆管实验中，圆管内径D为5.87 mm，圆管的加热段为60 D，壁厚为0.28 mm。方法验证中采用的实验工况如表1所示。验证工况包含了高热流密度和低热流密度。本文以氦氙混合气体为工作流体，其物性采用Tournier等人[16]提出的方法计算，图4给出了3.0 MPa下氦氙气体物性随温度的变化图。

Figure 3. Experimental section of Taylor’s round tube

图3. Taylor圆管实验段

Figure 4. Variation diagram of He-Xe physical properties with temperature

图4. He-Xe物性随温度变化图

图5给出了实验数据与计算数据的对比。从图中可以看到，文中使用的数值方法与实验值符合较好。在圆管计算中，计算值与实验值在低热流密度下误差在3%以内，在高热流密度下误差在5%以内。

2.4. 计算参数设置

在紧密栅通道的研究中，边界条件设置如表2所示，通道的入口和出口分别定义为速度入口(或者质

Table 1. Boundary condition settings of the simulation example

表1. 模拟算例的边界条件设置

Figure 5. Taylor’s experimental simulation results: (a) Wall temperature in working condition 1; (b) Wall temperature in working condition 2

图5. Taylor实验模拟结果：(a) 工况1壁面温度；(b) 工况2壁面温度

Table 2. Boundary condition settings

表2. 边界条件设置

量流量入口)和压力出口，这些边界条件设置参考了大功率高温气冷的边界条件并做出了一些的调整。通道的包壳内壁面采用均匀热流密度进行加热，外壁面被定义为无滑移边界条件并采用耦合模式，六角形管道壁面被定义为无滑移边界条件并施加绝热边界条件。流体为40 g/mol的He-Xe混合气体。模拟过程中压力–速度耦合采用SIMPLE格式，连续性、动量、能量等的收敛准则为小于10⁻⁶。为了保证数值模拟的收敛性，对出口处的平均截面压力进行了监测。

3. 紧密栅结构下的流动换热特性

在工况1的情况下，计算了He-Xe混合气体19棒模型(六边形绝热边界)中的流动换热。图6给出19棒模型距离入口0.5 m处的轴线速度、横向速度和出口的温度分布图。从图6中的云图可见中心区域和外围区域以及边缘区域子通道中速度温度场的区别。从图6(a)中可见，中心区域速度较低于边缘区域速度，这是因为在紧密栅结构下，在中心棒束狭小的通道之间流动存在较大局部阻力，而边缘区域通道间隙较大。从图6(b)中可见，随着轴向的不断发展，中心区域和边缘区域的温度差异逐渐变大，中心区域温度较高，边缘区域温度较低，这是因为边缘区域的冷却剂流速较大，能更快地带走热量。同时，棒束表面的温度存在不均匀的分布，燃料棒间隙之间的温度较高。从图6(c)可见，19棒通道中存在一定的横向流动以及交混，最大横向流速为0.486 m/s，横向流动从中心区域流向边缘区域，产生这种横向流动的主要原因是区域之间的流动差异引起的压差。从图6(c)可见，在入口段的横向流动较强，在充分发展后趋于一致，这是因为流体在入口段还处于发展阶段，还没有形成明显的流动边界层。

Figure 6. Velocity and temperature distribution diagrams: (a) Velocity distribution cloud diagram; (b) Temperature distribution cloud diagram; (c) Transverse velocity distribution cloud diagram

图6. 速度温度分布图：(a) 速度分布云图；(b) 温度分布云图；(c) 横向速度分布图

为量化不同位置子通道之间的流动换热差异，选取中心区域、外围区域以及边缘区域的子通道，给出了沿着y方向上的速度温度分布。图7为充分发展段(距离入口0.5 m)平面上不同区域子通道的流体温度速度分布。从图7(a)可见，不同子通道之间的速度差距不大，边缘区域的棒间隙速度大于中心区域，这是因为边缘区域子通道仅有三根燃料棒作为约束，而中心区域和外围均为4棒作为约束。从图7(b)可见，中心区域的棒间隙温度高于外围以及边缘区域，这是因为外围区域以及边缘区域的通道均与流速较大温度较低的边缘子通道相连。从壁面温度来看，中心区域和边缘区域相同，外围区域的壁面温度较低，猜测是因为外围区域的棒束在边通道的高流速冷却以及包壳的导热下，导致整体壁面温度较低。整体而言，中心区域的温度整体上要大于外围区域和边缘区域。

4. 不同栅径比结构下的流动换热特性

在紧密栅棒束通道的流动换热中，栅径比有着重要的影响，探究He-Xe混合气体在什么栅径比的紧密栅通道中具有较好的换热性能具有重要意义。本节对P/D = 1.1、1.15、1.2、1.25、1.3五种不同栅径比的流动换热进行了研究，这五种栅径比紧密栅通道仅改变其栅距，在相同的边界条件以及入口雷诺数(Re = 4.2 × 10⁴)工况下进行栅径比对He-Xe混合气体流动换热特性的影响的研究。比较了不同栅径比下的局部努塞尔数、范宁摩擦系数以及综合换热性能，它们的计算式如下：

$N{u}_{local}=\frac{q{d}_h}{\lambda \left(\overline{T_w\left(z\right)}-\overline{T_b\left(z\right)}\right)}$ (9)

$f=\frac{2\Delta p{L}_c}{\rho L{v}^2}$ (10)

$THRR=\frac{Nu}{N{u}_{ref}}/\frac{f}{f_{ref}}$ (11)

图8为相同雷诺数工况(工况2)下，不同栅径比下的局部努塞尔数、范宁摩擦系数随着轴向距离的变化图。从图8(a)中可见，Nu数随着轴向距离的变化趋势均为先快速减小后增加再缓慢减小，猜测这是

Figure 7. (a) Distribution diagram of fluid velocity in different subchannels; (b) Distribution diagram of fluid temperature in different subchannels

图7. (a) 不同子通道中的流体速度分布图；(b) 不同子通道中的流体温度分布图

Figure 8. Local Nusselt number (a) and friction coefficient along the path (b) at different pitch-diameter ratios

图8. 不同栅径比下的局部努塞尔数(a)和沿程摩擦系数(b)

因为在入口段流体还处于发展阶段，边界层较薄，热阻较小，因此对流换热较强；随着轴向流动的发展，边界层变厚，对流传热变弱；随着流动从层流发展为湍流，对流传热增强；在之后的充分发展区域，由于壁面温度不均匀分布的加剧，棒间隙之间的流动换热变差，对流传热变弱。随着边界层的发展，从图8(b)中可见，范宁摩擦系数沿程逐渐减小，这是因为在流动过程中，随着摩擦损失，流体的平均流速降低，摩擦系数减小。从图8(a)可见，在入口处，局部努塞尔数随着栅径比的增大而增大；在充分发展段时，局部努塞尔数随着栅径比的增大而先增大后减小，在P/D = 1.25时达到最大。从图8(b)可见，随着栅径比的增大，范宁摩擦系数也呈现先增大后减小，在P/D = 1.25时达到最大。图9(a)给出了平均努塞尔数和平均摩擦系数随栅径比的变化曲线。在五种栅径比中，P/D = 1.25的平均努塞尔数最大，其平均摩擦系数也最大。用性能评价标准对5种不同栅径比的紧密栅棒束通道中综合换热性能进行了评价，结果如图9(b)所示，图中可见随着栅径比的增大，紧密栅棒束通道的综合换热性能先增大后减小，5种栅径比中P/D = 1.15的综合换热性能最佳。

Figure 9. Average Nusselt number and average frictional factor (a) and comprehensive heat transfer performance (b) at different grid size ratios

图9. 不同栅径比结构的平均努塞尔数和平均摩擦阻力系数(a)，以及综合换热性能(b)

5. 总结

紧密栅气冷空间堆堆芯结构具有功率密度高、质量小、结构紧凑等优点，本文对40 g/mol的He-Xe混合气体在P/D = 1.2的紧密栅棒束通道中的流动传热特性进行了数值分析，对其温度场、速度场以及横向流动特性进行了定性和定量的描述。此外，比较了5种不同栅径的比紧密栅结构的流动换热特性，计算了在同一雷诺数工况下的局部努塞尔数和摩擦系数，得到了不同栅径比下的综合换热性能。

1) 构建了适用于He-Xe气体流动传热计算的CFD模型，其中湍流模型采用修正后的SST k-ω模型，物性计算采用Tournier方法，计算结果与实验结果误差不超过5%。

2) 得到了He-Xe在P/D = 1.2的19棒紧密栅结构下的流场、温度场分布以及局部努塞尔。结果显示，中心区域的子通道的速度分布相似，差异不大；中心区域的温度整体高于外围和边缘的温度。紧密栅机构下的Nu数随着轴向距离的变化趋势为先快速减小后增加再缓慢减小。

3) 比较了不同栅径比(P/D = 1.1~1.3)下，氦氙混合气体的全局努塞尔数和摩擦阻力系数的变化；用性能评价标准对5种不同栅径比的紧密栅棒束通道中的综合换热性能进行了评价。结果表明：综合性能评价标准随着栅径比的变化先增加后减小，换热性能最佳的栅径比为P/D = 1.15。

致 谢

感谢周源老师对文章写作的指导。

NOTES

^*通讯作者。

参考文献