温永江，王成龙，张大林，秋穗正，苏光辉，田文喜

(西安交通大学 核科学与技术学院，陕西 西安 710049)

核热推进反应堆具有超长续航、机动性强、隐蔽性高、突破能力强等巨大优势，2018年俄罗斯“海燕号”的公布更是将世界的目光聚焦到核动力领域[1]。但核热推进反应堆会在运行后产生大量的放射性物质，同时出于实用性考虑，核热推进反应堆没有完善的辐射屏蔽装置[2]，会对外界环境造成大量的辐射污染[3]。

目前辐射场计算主要有以下几种方法：离散纵标法、蒙特卡罗方法和点核积分法[4]。这3种计算方法都有各自的优点和缺点。离散纵标法可将三维转化为一维或二维问题，快速求解几何较简单的问题[5]；蒙特卡罗方法则采用概率论的方法，计算结果较准确同时需要的计算时间也较长[6]；点核积分算法则采用半经验公式，计算速度快[7]。目前研究者普遍关注将3种方法耦合来加速计算。郭雨菲等[8]开展了将蒙特卡罗方法与点核积分算法耦合来进行快速计算的研究，降低结果误差同时大幅提高计算速度，宾西法尼亚大学开发的A3MCNP程序使用离散纵标加速MC法，可初步解决深穿透问题。目前国内外针对点核积分算法展开了大量研究，研发出QAD-CG、QAD-CGA、GSCM、PKI等程序[9]，点核积分算法已经成为辐射屏蔽计算的基本方法之一。目前世界上已经公布的核动力导弹型号有两种：俄罗斯的“海燕”号[10]和美国的“冥王星”号[11]。其中“冥王星”资料更公开透明，作为本文研究对象。1957年，美国开始实施“冥王星”计划，最终研发出核热推进导弹“SLAM”，先后提出两代核热推进反应堆“TORY Ⅱ-A”和“TORY Ⅱ-C”。其中二代堆型“TORY Ⅱ-C”功率更大，设计更成熟，本文选择“TORY Ⅱ-C”反应堆为计算对象。本文主要针对TORY Ⅱ-C运行后产生的放射性元素造成的辐射进行研究，目的是为反应堆运行结束后回收时提供辐射屏蔽建议。

1 几何建模和计算方法

1.1 几何建模

美国Lawrence辐射实验室(LRL)20世纪提出的TORY Ⅱ-C靠空气冷却堆芯[12]，UO2-ZrO2-Y2O3颗粒弥散在BeO基体中[13]，直径通常为0.5～1 μm，表1列出了堆芯的主要参数[14]。图1和表2分别给出了该堆的剖面图及燃料排布，径向分18层，轴向分6层。在燃料区域前后以及侧面都存在反射层。前反射层为25 cm厚BeO，后反射层为6 cm厚BeO，侧面反射层为5 cm厚BeO加外层2.5 cm厚Ni。堆芯采用六棱柱燃料组件，图2示出了模拟的六棱柱燃料组件和堆芯的建模图，其中每种颜色都代表1种UO2质量分数不同的燃料。

表1 TORY Ⅱ-C主要参数Table 1 Main parameters of TORY Ⅱ-C

图1 堆芯建模图Fig.1 Modeling diagram of core

表2 燃料分布Table 2 Fuel distribution

图2 燃料组件和堆芯建模图Fig.2 Modeling diagrams of single fuel assembly and core

1.2 计算方法

同其他反应堆相同，在不考虑核热推进反应堆排出气体的情况下，TORY Ⅱ-C对外的辐射源项有3个：1) 链式裂变反应导致的直接辐射；2) 反应堆构件活化后的放射性；3) 裂变反应后放射性产物的辐射。本文假设的场景为核热推进反应堆经过运行后回收处理的过程，此时反应堆已经停止运行，链式裂变反应的直接辐射不再考虑。与核电厂反应堆不同，核热推进反应堆具有能量密度高、体积小、重量轻的特点，为实现小型化设计，提高核反应性，核反应堆必然采用高浓缩铀材料，放射性运行产物将是停堆贮存期间的主要放射性源项，本文针对TORY Ⅱ-C运行后产生的放射性产物造成的辐射进行研究，分别进行了燃耗计算、能谱计算以及辐射剂量计算。本文在进行燃耗计算时，将蒙特卡罗算法和燃耗方程进行耦合，分别求解中子输运方程与燃耗方程。

中子输运方程：

∬Σs(r,E′)f(r,E′→E,

∬vΣf(r,E′)Φ(r,Ω′,E′)dΩ′dE′

(1)

燃耗方程：

(2)

之后根据燃耗计算结果计算光子能谱和中子能谱，计算的伽马射线源强度和能谱包括X射线、伽马射线、韧致辐射、自发裂变伽马射线和伴随(α，n)反应的伽马射线产生的光子，中子源强度和能谱包括自发裂变(α,n)反应和(β)反应。

计算辐射剂量时采用点核积分算法通过求解方程(3)，计算反应堆外环境剂量率：

(3)

计算辐射剂量时将堆芯结构进行了简化，堆芯中燃料与空气流道中的空气均匀化，视作一种物质。TORY Ⅱ-C燃料管孔隙率约53%，再加上控制棒通道等，燃料区域孔隙率约能达到57%，将燃料密度倍增系数设置为0.43即可，简化后的模型如图3所示，包括均匀化的燃料区域、前反射层、后反射层以及侧反射层。通过这种方式来构建几何模型之后，反应堆结构变为中心对称，二维辐射场可视化可展示剂量当量率的分布情况。

图3 简化模型Fig.3 Simplified model

2 结果与讨论

2.1 源项计算

TORY Ⅱ-C反应堆最大功率为500 MW，设计运行时间为10 h，本文将计算工况定为反应堆满功率运行10 h。计算结果包括全堆运行产物的核素种类和放射性活度，以及各燃料区域的核素种类和放射性活度，获得更详细的源项分布结果有利于之后进行精细的辐射屏蔽计算。

通过计算得到，235U从0.081 75 g/cm3降到0.081 29 g/cm3，共消耗0.56%。放射性活度随时间的增加逐渐变大，最后达到1.298×1021Bq，如图4所示。为方便后续计算，通过对锕系产物和裂变产物放射活性进行分析，共筛选出106种放射性较高的核素进行后续能谱和辐射屏蔽计算，其放射性活度达到1.09×1021Bq，占到锕系产物和裂变产物放射性活度的约90%，保证了一定的精确性。计算时没有设置衰变时间，106种核素中存在很多衰变时间很短的核素，若添加几个小时或几个月的衰变时间，结果将很不同，同时放射性活度也将快速下降。和裂变产物相比，锕系产物放射性活度很小，占比不足0.1%，在后续辐射计算中106种核素只有裂变产物。通过与功率分布比对，发现裂变产物分布与之匹配。TORY Ⅱ-C燃料排布很细，分布合理，周边也存在反射层，功率分布与裂变产物的分布都比较均匀。

图4 反应堆放射性强度随时间的变化Fig.4 Temporal variation of radioactivity

针对放射性活度较高的106种核素进行计算，获得其中子能谱与光子能谱，为后续辐射剂量计算提供源项输入。中子粒子数量总共为2.316×1018s-1，光子粒子数量总共为2.482×1021s-1，光子数量是中子的1 000多倍，中子辐射占比非常小，以伽马辐射为主，在后续辐射计算中只进行伽马辐射计算。

2.2 辐射计算

输入辐射源分布参数时，将反应堆沿轴向分为6层，层与层之间能谱和源强不同，层内通过设置权重调整径向的源分布，之后分别计算出每层的辐射强度，再将每层计算结果相加即可得到最终的剂量率，这样可以模拟出放射性物质在堆内不均匀的分布特征。

图5为反应堆周围剂量率的分布情况，白色区域为反应堆本体。反应堆尾部，坐标(0，-6)，为外环境剂量率最高的位置，剂量当量率达到67.185 Sv/h。相比较而言，反应堆正前方，坐标(0，154)剂量率为15.152 Sv/h，比后方区域的剂量率低，其主要原因是反应堆前后反射层结构不同，反应堆前反射层为25 cm厚BeO(孔隙率57%)，而后方只有6 cm厚BeO(孔隙率57%)。反应堆侧面剂量率普遍更低，坐标(68，90)处剂量率只有33 mSv/h，坐标(68，-6)处剂量率只有8.38 mSv/h，坐标(68，154)处为2.85 mSv/h，因为侧面反射层为5 cm厚BeO(孔隙率10%)加外层2.5 cm厚Ni(孔隙率10%)，Ni的辐射屏蔽效果极佳。若将Ni替换成等厚度的BeO，3个位置的剂量率分别为28.987 Sv/h、10.493 Sv/h和3 Sv/h，也会达到很强的辐射水平。结合表3所示某快堆辐射屏蔽设计的相应标准[15]，可看出运行后的反应堆对环境的辐射是很高的，远超0.5 mSv/h。

图5 反应堆周围辐射场Fig.5 Radiation field around reactor

表3 某快堆辐射屏蔽设计标准Table 3 Design criteria for radiation shielding of some civil fast reactor

图6示出了将放射性物质视作非均匀分布与视作均匀分布的中轴线位置计算结果对比。比较结果来看，将源视作均匀分布计算时得出的结果与非均匀分布的计算结果相差不大。视作均匀分布时，前侧的计算结果普遍高10%左右，尾部普遍低3%。原因是后侧放射性物质比平均水平高，而前侧则较低。图6也显示出随着计算点与堆距离的增加，结果偏差的变化。计算点与反应堆越远，结果偏差越小。综合来看，若偏差在可接受范围内，可为了计算效率将放射性物质看作均匀分布或适当减少堆的计算区域数量，采取的计算区域数量可进行比对计算选择合适的数值。

图6 放射性核素均匀分布与非均匀分布的剂量率对比Fig.6 Comparison of dose rates between uniform and non-uniform distributions of radionuclides

图7示出了反应堆周围小范围的剂量率分布情况，白色区域为反应堆本体。图8为表面剂量率随半径的变化，可看到位置越靠近堆的中心轴，该点的剂量率越高。随着半径的增加，剂量率减小的趋势越来越大。在半径65～67.5 cm处剂量率发生突变，这与侧面反射层出色的伽马辐射屏蔽能力有关。

图7 反应堆周围小范围辐射场Fig.7 Radiation field around reactor for small scale

图8 表面剂量率随半径的变化Fig.8 Change trend of surface dose rate with radius

图9为侧反射层外表面剂量率变化。可看出，整体呈现中间高，两边低的趋势。轴向坐标10～110的区域内，剂量率达到了很高的水平，燃料区域范围为0～131。两侧剂量率低主要有两个原因：1) 因为燃料区域边缘辐射物质放射性活度低；2) -6～0、131～154范围内没有燃料和放射性物质，只有反射层。侧面屏蔽防护应以燃料所在区域的中心部分为重点。

图9 侧反射层外表面剂量率变化Fig.9 Change of dose rate on outer surface of side reflector

3 结论

本文针对美国TORY Ⅱ-C回收后运行产生放射性裂变产物造成的辐射展开研究，进行了源项计算以及辐射计算，最终得到该种反应堆满功率运行10 h后的辐射场，主要结论如下。

1) 10 h满功率运行共消耗0.56%的235U。其中裂变产物的放射性活度远超锕系产物，达到1.298×1021Bq。裂变产物释放的辐射粒子中光子占到了99.9%以上，反应堆以伽马辐射为主。

2) 运行后反应堆辐射水平很高，反应堆屏蔽层外表面最大辐射剂量达到67.185 Sv/h，位于反应堆后反射层外表面中心，此值远超人体可以接受的水平。侧面因为存在2.5 cm厚的Ni，剂量当量率大幅减小，侧反射层外表面最高值只有33 mSv/h。比较而言，尾部辐射最高，头部次之，侧面最小。同时反应堆中心轴位置普遍比周边区域剂量率更高。

3) 将辐射产物视作均匀分布与分区域非均匀分布的计算得到的结果相差不大，在较近距离内差值最大为12%，且随着与反应堆距离的增加，差值还会继续减小。

综合上述计算分析，本文提出以下建议。

1) 核热推进反应堆没有完善的辐射屏蔽装置，必须在回收存放时增设一定厚度的屏蔽设施降低其放射性影响，其中反应堆中轴线及尾部区域是屏蔽的重点位置。反应堆以伽马辐射为主，应选用原子序数高、密度高的材料作为屏蔽材料，如铅、铁和混凝土等。

2) 计算剂量当量率时，为提高计算效率，可将源的分布做简化处理甚至可视作均匀分布。这种方法在计算距离反应堆较远位置的剂量率时十分有效。