空间气冷反应堆堆芯流动换热数值仿真研究

孟 涛，赵富龙，程 坤，曾 陈，谭思超

(哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001)

近年来，诸如空间反应堆等特种反应堆开始受到越来越多的关注。此类反应堆因功率较小、结构相对简单，利于反应堆的小型化，因而在对空间体积和重量要求较为严格的空间环境、船舶环境中具有更为广阔的应用前景。其中直接布雷顿循环气冷反应堆因结构紧凑、体积重量较小、功率适用范围大等特点而多应用于大功率空间系统[1]。

美俄均提出过多个空间气冷反应堆研究计划[2-4]。2003年3月，美国航空航天局签署了普罗米修斯核能系统研究计划[1]，该计划作为木星冰卫星探测器计划的后续计划，目标在于研制一款安全、可靠、能执行多种外太阳系探索任务的深空探测器。该探测器以核电推进为动力，可携带任务载荷不小于1.5 t，并要求任务相关技术可应用到月球与火星表面供电任务以及空间货物运输任务上。虽然该任务开展仅两年后便被终止，但该任务对空间反应堆相关技术进行了广泛的调研，针对空间用反应堆堆芯类型、热电转换方式、反应堆结构材料等，给出了全面、详细的技术调研报告，能为相关研究提供很好的帮助。另一方面俄罗斯也提出了兆瓦级空间核反应堆计划[5]，该计划旨在开发兆瓦级功率空间供电与推进反应堆系统，能为空间探索及输送任务提供动力。

为揭示空间气冷反应堆堆芯流动换热机理，本文基于美国普罗米修斯计划的公开信息[1,5]，建立相应的空间反应堆堆芯CFD模型，针对空间反应堆堆芯流动与传热特性开展计算分析，评估现有设计中仍待优化之处，并提出相关的优化建议，为未来空间气冷反应堆设计提供参考。

1 计算几何模型

普罗米修斯计划提出了多个气冷堆方案，包括基准方案、低功率方案、开放流通区域方案等，但计划终止时并未给出最终设计方案[6]。为分析空间反应堆的流动换热特性，本文参照基准设计方案进行数值仿真计算。反应堆堆芯内燃料棒以正六边形排布，由内至外共排布有7圈燃料棒，最外层1圈燃料棒在正六边形的6个角各缺少1根燃料棒，燃料排布如图1所示。中心区域不排布燃料棒，仅设置1根控制棒，且控制棒通道与冷却剂通道隔离，堆芯几何模型如图2所示。

图2 普罗米修斯计划空间反应堆结构

反应堆由内、外两层压力容器构成，内外两层压力容器通过开有孔洞的栅板连接。燃料棒结构及冷却剂流动方式如图3所示。每根燃料棒具有独立的环形冷却剂通道，燃料棒通过入口、出口两个端部与环形冷却剂通道及堆芯上下支撑板固定。堆芯主要参数信息列于表1。

图3 燃料棒结构

参数数值燃料棒外径,mm19.65冷却剂通道厚度,mm2.16燃料棒长度,m1.118燃料棒芯块直径,mm18.19燃料活性区长度,m0.608燃料棒数量288堆芯压力容器外径,mm618.1压力容器内环隙厚度,m0.3冷却剂质量流量,kg/s6.16堆芯入口温度,K911堆芯出口温度,K1150反应堆热功率,kW1002

由于文献中的堆芯几何参数信息有限，因此本文使用的部分堆芯结构参数为推测数据，包括内压力容器内径、冷却剂主管道内径和燃料棒出口端结构尺寸。数据的推测依据为普罗米修斯计划反应堆结构的几何比例。推测数据对堆芯内部流动与换热影响较小，主要影响进出口处的流速及流场特征。

2 计算模型及条件设置

2.1 边界条件设置

由于该反应堆的冷却剂通道为独立通道，故选取单根燃料棒及其冷却剂通道作为单通道模型进行计算。首先建立了单通道模型进行仿真，以进行网格无关性验证。入口条件均设置为质量流量入口，经过计算确定单根燃料棒通道平均冷却剂流量为0.021 389 kg/s，相应的1/6堆芯冷却剂流量为1.026 67 kg/s，入口温度为911 K，出口条件设置为压力出口，压强为1 898 kPa。单通道燃料棒加热方式为余弦加热，积分加热功率为3.479 kW。1/6堆芯的加热方式采用真实堆芯功率分布，堆芯出口平均温度控制在1 150 K。堆芯功率分布通过使用SuperMC蒙特卡罗软件[7-8]计算得到，如图4所示。在热工水力计算中使用线性表格插值的方式调用功率分布。计算中使用He-Xe混合气体作为冷却剂[9-11]，混合摩尔质量参照普罗米修斯方案选定为31.5 g/mol。

图4 反应堆堆芯功率分布

2.2 网格模型

多面体网格与其他非结构网格相比具有最少的网格数量、较少的计算时间和较高的计算精度，因此本文采用多面体网格，网格划分采用自适应加密方法，可根据结构特点自动加密网格，基准网格尺寸为6 mm，最小网格尺寸为基准网格尺寸的25%，近壁面区采用边界层网格，边界层网格层数为5，总厚度为1.2 mm，最终划分的局部网格结构如图5所示。

2.3 物理模型

由于目前关于He-Xe混合气体的相关实验公开文献较少[12-16]，且实验精度难以保证，因此，本文使用数值计算方法模拟空间气冷反应堆堆芯冷却剂流动情况，并参照文献[15]中的对流换热实验数据，进而获得堆芯流动换热特性，为空间反应堆的优化设计提供指导。本文使用适用范围较广的可实现k-ε模型进行数值计算。可实现k-ε模型是两方程雷诺平均湍流模型，由Shih等[17]基于标准k-ε模型推导得出。由于该模型对湍流黏性系数的计算采用了时均应变率、旋度及旋转坐标系角速度的函数，因而相比于标准k-ε模型通常具有更好的计算精度和网格适应性。

图5 网格结构

3 计算结果及分析

3.1 网格无关性验证

图6 参数随网格数量的变化

为确保计算结果的精度和可靠性，首先需进行网格无关性验证。图6示出了利用单通道模型计算得到的网格敏感性分析结果。结果表明，当网格数量大于30万时，冷却剂出口最高温度与平均速度的变化幅度开始显著变小，相对误差在3%以内。因此，对于单通道，网格数量应控制在30万以上。而对于1/6堆芯，由于其包含54根燃料通道、上下腔室以及下降段，因此本文选择的网格数量为2 558万，网格质量大于0.9。

3.2 计算模型验证

为确保计算模型的可靠性，需进行计算模型验证。由于使用He-Xe混合气体作为冷却剂的流动换热特性及仿真属于细分研究领域，且实验成本高昂，目前为止公开发表的实验文献及仿真研究极其有限[12-16]。本文参照实验参数相对全面的文献[15]进行计算模型验证。计算使用可实现k-ε模型，网格采用多面体网格，网格数量为52.5万，计算边界条件列于表2，计算结果示于图7。由计算结果发现，冷却剂中心温度与文献[15]中实验数据符合较好，最大温度误差小于5%。但加热管壁面温度相差较多。通过在流固交界面增加接触热阻0.005 m2·K/W后，加热管壁面温度在整体趋势上与文献[15]中的实验数据能较好地符合。由此说明，使用可实现k-ε模型可较好地模拟He-Xe混合气体管内流动的流体特性，但仿真得到的壁面温度会偏低，进而导致计算得到的努塞尔数Nu或对流换热系数h偏低。增加接触热阻后，可对壁面温度进行较好的仿真模拟。由于本文进行的1/6堆芯计算，未考虑流固耦合，故未在计算中增加接触热阻。

表2 计算模型验证算例参数

图7 仿真验证算例与文献实验数据对比

3.3 堆芯流场分析

为简化计算，忽略燃料棒、控制棒、格架等结构部件，仅针对流体区域进行数值计算。通过计算得到的冷却剂在入口段和出口段的速度云图如图8所示。由图8a可知，在入口段，流体进入压力容器间环隙后直接撞击在内层压力容器壁面，形成分流，并在撞击点附近形成低流速区。向上方弧形段流动的冷却剂流速快速降低，向下方流动的冷却剂通过开孔栅板时由于流通面积的减少而出现加速流动。通过开孔栅板后流速逐渐降低，流动到下方开孔栅板时又出现加速流动，之后流速再逐渐降低，并开始向堆芯中心侧偏转，如图8b所示。在偏转流动过程中，靠近弧形壁面的一部分入流冷却剂继续沿壁面向下流动，直至流动到堆芯纵向中轴线位置处，在受到中心控制棒通道壁面阻挡后改变流动方向，沿控制棒通道壁面竖直向上运动，与入流冷却剂混合后共同流入燃料棒冷却剂通道。流出燃料棒通道后，不同冷却剂通道的流体汇聚后通过出口管流出反应堆。

图8 冷却剂入口段(a)和出口段(b)速度场

图9 出口管横截面速度场

为了解冷却剂在出口位置处的速度分布情况，图9示出了出口管横截面的速度云图。由图9可知，在冷却剂流出管横截面上，速度分布不均匀，上方速度较高，下方速度较低。

图10示出了不同高度位置处的堆芯横截面速度分布结果。在入口管中间位置(图10a)处，堆芯速度场沿出口管中间轴线呈对称分布。在圆周方向中间位置处，冷却剂流速较高，而靠近内压力容器壁面处的外围冷却剂流速较低。冷却剂流出燃料棒区域后快速进入出口管，并未在上封头内大范围流动，因而使得流速分布不均匀，产生了局部高流速区和低流速区。在入口段下方约8 cm位置(图10b)处，环状下降段内冷却剂流速分布并不均匀，入口管内的冷却剂进入后受压力容器壁面的影响形成分流，导致在径向方向上冷却剂速度分布不均匀，进而影响沿环状下降段流动的冷却剂的速度分布。

为了解堆芯通道中的温度分布情况，图11示出了冷却剂的温度分布云图。由图11可知，当冷却剂进入堆芯后，由于径向加热的不均匀性，堆芯径向中间位置处冷却剂具有最高温度，而半径较小和较大位置处的冷却剂温度都相对较低。当冷却剂流入出口端下封头时，冷却剂并未发生明显搅混，进而造成该位置处冷却剂温度分布不均匀。

图10 入口管段中间位置(a)及下方8 cm处(b)冷却剂速度场分布

图11 堆芯温度场分布

此外，为了解冷却剂在堆芯通道中流动过程中流动阻力的变化情况，为结构优化设计提供参考，图12示出了堆芯通道内的压强分布云图。由图12可知，冷却剂在流动过程中的压降主要发生在3个位置，1个是环状下降段的2处支撑处，另1个则是在堆芯燃料区域的冷却剂通道中。

图12 堆芯压强分布

堆芯入口处，由于入口管道几乎垂直于内压力容器，因此入流冷却剂直接撞击内压力容器，产生较大能量损失。从优化的角度考虑，可改变入口管道方向，使其向上方倾斜，与压力容器下降环段呈钝角，以避免入流冷却剂直接撞击内压力容器，进而减少不必要的能量损失。

另一方面，由图8a可知，冷却剂流入压力容器后，一部分向上方的弧段环隙流动。但弧段环隙内大部分流体流动较为缓慢。因此，可考虑对上部弧形封头进行结构优化设计，减轻系统质量。

此外，在冷却剂出口段，冷却剂搅混程度较轻，导致出口段管道、弧形封头内存在较大的温差与速度差，不利于平衡封头结构和出口管道的内部热应力，既有可能降低结构寿命，同时还会影响气轮机工作性能。因此，从优化的角度考虑，可考虑设置加强搅混结构或优化出口管排布位置，来达到降低冷却剂出口温度分布不均匀的目的。具体优化方式仍有待进一步计算分析。

图13 不同冷却剂通道内温度分布

为评价不同冷却剂通道内流体的流动特征，图13示出了不同冷却剂通道内冷却剂温度随堆芯高度的变化，图中标号如图14所示。由图13可知，不同通道的温度变化均呈现相同的升温趋势，升温速率均是先小后大，然后再逐渐减小。越靠近半径中间位置，加热功率越大，通道出口温度越高。计算结果显示，不同通道出口温度最大温差可达76 K。因此，当堆芯出口处搅混不够充分时，堆芯出口管内冷却剂将产生较大的温度不均匀分布，进而影响气轮机设备的正常运行。

图14 冷却剂通道示意图

3.4 不确定性分析

本文参照文献[18]中的不确定分析方法，针对3种不同网格情况进行仿真结果的不确定性分析，网格数量分别为2 558万(网格1)、4 284万(网格2)和4 733万(网格3)。计算得到的冷却剂出口温度、平均流速和压降列于表3。根据计算结果得到的网格收敛系数(GCI)列于表4，其中网格1/2表示利用网格1与网格2计算结果得到网格收敛系数，依次类推。将最大网格收敛系数除以2即可得到数值不确定度。由此可知，仿真计算的冷却剂出口温度不确定度为0.908 5 K，冷却剂出口流速不确定度为0.182 m/s，反应堆压降不确定度为1 419 Pa。

表3 不同网格数量下的计算结果

表4 网格收敛系数

4 结论

本文针对美国普罗米修斯计划提出的反应堆方案建立了三维模型，并使用可实现k-ε模型进行了1/6堆芯的流动换热计算，堆芯加热功率使用堆芯稳态物理计算得到的功率分布曲线，通过计算得到了堆芯温度场、速度场和压力场的分布情况，并分析了空间气冷反应堆堆芯的流动换热特性，得到以下主要结论。

1) 在入口段，流体直接冲击内压力容器壁面，形成分流，并在撞击点附近形成低流速区。从优化的角度考虑，可改变入口管道方向，使其向上方倾斜，与压力容器下降环段呈钝角，以避免入流冷却剂直接撞击内压力容器，进而减少不必要的能量损失。

2) 在堆芯下降段及燃料区域，冷却剂速度场沿出口管中间轴线呈对称分布。上腔室内圆周方向中间位置处冷却剂流速较高，而两侧靠近内压力容器壁面处的冷却剂流速较低，环状下降段内冷却剂流速分布并不均匀，入口管内的冷却剂进入后受压力容器壁面的影响形成分流，导致在径向方向上冷却剂速度分布不均匀，进而影响沿环状下降段流动的冷却剂的速度分布。

3) 堆芯出口管截面上温度及速度分布不均匀，温差可达76 K，从优化的角度考虑，可设置加强搅混结构或优化出口管排布位置，降低冷却剂对气轮机等设备的影响。