0.   引　言

微型反应堆具有寿命长、供电稳定、对外界的依赖性小等优点，可为偏远基地、空间探索等远程环境的应用提供可靠的移动电源，是未来深海、深空和陆上应急供电等领域特种能源的最优选择，有重大的潜在利用价值，有可能成为改变未来核动力格局的颠覆性技术之一。微型反应堆在动力学冲击发生的情况下会发生剧烈形变，从而引发反应性变化，如何准确快速计算撞击引入的反应性变化成为了微型反应堆设计和安全分析中亟待解决的关键问题，比如空间核反应堆电源的使用中就要求在达到工作轨道前发生一切可能事件时均不能进入临界状态^[1]，最常用的潜在反应性增加场景是反应堆在发射台附近的混凝土地面上发生撞击。空间核反应堆在撞击的情况下会发生剧烈形变，从而引发反应性变化。

目前国内还没有研究冲击对反应堆物理临界计算的研究，研究人员大部分仅对掉落到水或湿沙的事故环境影响进行了评估^[2-4]，没有对撞击事故可能造成的形变进行模拟分析，可能原因是撞击导致的核反应堆或燃料棒束的反应性变化过程很难精确建模，尤其是涉及到大变形条件下的多物理场耦合计算。在国外，橡树岭国家实验室采用桑迪亚国家实验室研发的CTH程序对空间核反应堆进行撞击地面的模拟^[5]，但是该方法目前只适用于二维堆芯模型，无法进行复杂三维建模。桑迪亚国家实验室采用有限元程序PRESTO和DAG-MCNP程序耦合的方法对空间堆的地面撞击事件进行安全分析^[6]，但该方法流程繁琐，DAG-MCNP程序在模拟反应堆形状畸变时，存在表面几何缺陷导致的粒子丢失、追踪错误和计算效率低等问题，一般需要特殊处理。

随着近些年连续介质力学分析程序（如ABAQUS）和并行计算机的发展，目前已经可以准确模拟反应堆高速碰撞下的大变形。在中子学计算方面，2012年蒙特卡罗程序MCNP开始支持非结构网格（UM）输运功能，经过后期不断发展，最新版MCNP的UM几何已经涵盖了组合实体（CSG）几何所有的物理输运问题，其主要将UM视为网格空间，将其添加到传统的MCNP的CSG几何中，在网格之外，以传统CSG几何对模型建模，网格内采用ABAQUS建模^[7-8]。得益于UM模拟的强大几何处理能力，国内外多家研究机构都开展了UM蒙特卡罗的相关研究工作，如清华大学的DAG-RMC、美国威斯康星大学及麻省理工学院的DAG-OpenMC和芬兰VTT技术研究中心的Serpent等。上述程序与MCNP相比都是表面UM，即把几何实体的表面进行面片化，对实体内部不划分网格^[9]。内部UM是把几何实体内部进行划分。与表面UM相比，内部UM在多物理耦合方面更有优势，可以与有限元程序采用相同的网格，可以支持基于网格的非均匀密度建模，更加符合真实情况^[10]，能准确地处理微型反应堆动力学冲击下的临界安全分析和新质堆研发及复杂几何辐射场生成问题^[11-14]。

基于上述研究，提出了一种基于非均匀密度的UM蒙特卡罗的方法来处理微型反应堆动力学冲击对反应性的影响，对代表性的反应堆装置进行初步计算分析，研究空间核反应堆在大变形情况下重返临界的可能性，为微型反应堆的应用安全评估奠定重要的理论基础。

1.   计算方法

为实现基于UM蒙特卡罗方法的空间堆动力学冲击引起的反应性变化的计算，本文基于Python脚本的方式建立空间堆动力学冲击反应性计算平台。在该平台中，基于ABAQUS程序完成瞬态下反应堆高速碰撞的模拟计算，再利用MCNP和DAGMC（本文选择DAG-OpenMC）计算变形堆芯的临界参数。空间堆碰撞引起的临界计算是通过每个时间步下模块间的调用实现的，图1为该耦合方法的计算流程。因为碰撞的过程很短暂，基本在几毫秒以内，从力学撞击导致的变形场数据传给临界计算过程中忽略了热工反馈等作用的影响，耦合过程目前仅仅考虑了中子输运和力学计算2个物理场，是单向耦合问题。当计算开始时，首先计算碰撞前堆芯有效增殖系数k_eff，然后基于特定的碰撞模型和接触边界条件，利用ABAQUS进行堆芯高速碰撞计算，再根据ABAQUS提供的每个时间步下变形后的网格更新MCNP和DAGMC的UM的密度，再进行输运计算。MCNP可以精细地考虑到同一材料不同网格内的非均匀密度变化，而DAGMC仅能考虑同一材料区均匀密度的变化。遍历完所有时间步长下的临界计算即可得到反应堆碰撞过程中k_eff、注量率分布等中子物理参数随时间的变化规律。

图  1  MCNP、DAGMC与ABAQUS耦合流程图

Figure  1.  Flow Chart of the Coupling of MCNP, DAGMC and ABAQUS

下载: 全尺寸图片 幻灯片

2.   计算模型

参考文献[12]，建立NaK冷却空间堆的堆芯模型，堆芯功率为25 kW，堆芯由85根UO₂燃料棒组成，²³⁵U的富集度达到93%，包壳是SS316不锈钢，采用NaK冷却剂，径向反射层选用Be，6个可旋转的控制鼓嵌入堆芯外的Be反射层中，控制鼓上有月牙形的B₄C吸收体，轴向反射层位于燃料棒的顶部和底部，选用BeO。堆芯最外层空白区域是一个水箱，装满硼化水，用于反应堆的径向屏蔽。堆芯的相关参数如表1所示。

表  1  NaK反应堆堆芯参数

Table  1.  Parameters of NaK Reactor Core

参数名	参数值		参数名	参数值
燃料棒数/根	85		燃料UO2密度/(g·cm−3)	9.873
燃料棒长度/cm	38		轴向反射层BeO密度/(g·cm−3)	2.85
燃料棒间距/cm	2.372		径向反射层Be密度/(g·cm−3)	1.85
燃料棒直径/cm	2		包壳SS316密度/(g·cm−3)	7.91
包壳厚度/cm	0.1		吸收体B4C密度/(g·cm−3)	2.52
反射层厚度/cm	15.3		围板SS316密度/(g·cm−3)	7.91
控制鼓直径/cm	15.7		压力容器SS316密度/(g·cm−3)	7.91
控制元件B4C厚度/cm	1.5		水箱容器SS316密度/(g·cm−3)	7.91

下载: 导出CSV 

| 显示表格

模型共分为7个部分，分别是堆芯组件、围板、压力容器、反射层、控制鼓、水箱和混凝土地面。碰撞计算中设置材料属性，包括材料密度、弹性、塑性。考虑到六面体网格可以提高计算结果的准确性和减少计算时间，划分网格时采用六面体网格对堆芯进行划分，其中单根燃料棒的网格数为4860，模型的总网格数接近140万。采用MCNP和DAG-OpenMC对反应堆进行建模计算，为了方便有限元分析，将包壳厚度设置为0.1 cm。建模时忽略了外围水箱容器厚度，只留下了硼化水的部分。

图2和图3分别展示采用CSG几何建模和采用ABAQUS的UM几何建模的NaK冷却空间堆几何图。采用ABAQUS模拟空间堆垂直和45°倾角跌落的情况，撞击的时间分别为0.0015、0.002 s，设置接触约束，将燃料（UO₂）、径向反射层（BeO）、包壳绑定在一起，混凝土靶标模拟过程中保持固定，设置空间堆的下降速度为100 m/s，核燃料UO₂的碰撞模型选择可压缩泡沫模型，该模型能够较好地处理核燃料撞击混凝土的问题，避免了对核燃料本身复杂物性和本构模型的依赖，与PRESTO采用的土壤泡沫模型类似，其他材料选择文献[12]中的本构参数，临界计算中考虑了混凝土的反射作用，混凝土采用了ABAQUS特有的CDP模型，应力-应变关系选用GB 50010-2010《混凝土结构设计规范》给出的混凝土本构关系，与文献[12]中PRESTO采用的K&C混凝土模型有一定差异，计算中采用ENDF/B VIII.0核数据库，投入粒子数10⁴个，考虑550代，忽略前50代，采用DAG-OpenMC的时候并没有选取与文献[12]一致的死单元删除方法，而是采用了允许一定比例粒子丢失的选项。

图  2  NaK反应堆的CSG建模

Figure  2.  CSG Model of NaK Reactor

下载: 全尺寸图片 幻灯片

图  3  NaK反应堆的UM建模

Figure  3.  UM Model of NaK Reactor

下载: 全尺寸图片 幻灯片

3.   计算结果分析与讨论

3.1   UM几何与CSG几何结果对比

表2展示了在控制鼓B₄C吸收体全部面向堆芯的情况下，MCNP采用CSG几何建模和UM几何建模以及DAG-OpenMC特征值的计算结果。可以看出，MCNP内部UM和DAG-OpenMC表面UM与MCNP的CSG几何下的特征值结果符合较好，MCNP-UM的偏差为129pcm（1pcm=10⁻⁵），DAG-OpenMC的偏差为101pcm，从而证实了本文UM蒙特卡罗计算的正确性。

表  2  UM与CSG几何下冷态临界k_eff的对比

Table  2.  Comparison of Cold Critical k_eff between UM and CSG　　　

堆芯状态	建模方式	keff	与CSG偏差
控制鼓面向堆芯	MCNP-CSG	0.87310
	MCNP-UM	0.87439	129pcm
	DAG-OpenMC	0.87209	101pcm

下载: 导出CSV 

| 显示表格

3.2   垂直跌落分析

图4为空间堆垂直跌落情况下ABAQUS模拟的最终位移结果。可以看出，应力分布存在明显的中心挤压现象，空间堆外层在与地面碰撞后中心出现塌陷现象。图5和图6展示了反应堆跌落前后几何和归一化注量率分布的变化情况，碰撞后注量率的最大值发生了70%以上的骤增，并且峰值位置下移。图7首先展示了垂直跌落特征值随时间的变化情况，可以看出MCNP-UM和DAG-OpenMC的k_eff计算值与文献[12]中SNL结果符合很好，差异主要来自本文的ABAQUS与SNL文献中PRESTO程序的不同，垂直撞击引起的k_eff增加约为8%，k_eff由碰撞前的0.87变为碰撞后的0.95。

图  4  垂直跌落位移云图

Figure  4.  Displacement Distribution of Vertical Impact

下载: 全尺寸图片 幻灯片

图  5  垂直跌落前后反应堆几何结构

Figure  5.  Reactor Geometry Structure before and after Vertical Impact

下载: 全尺寸图片 幻灯片

图  6  垂直跌落前后的中子注量率归一化分布

Figure  6.  Normalized Neutron Flux Rate Distribution before and after Vertical Impact

下载: 全尺寸图片 幻灯片

图  7  垂直跌落特征值和燃料区均匀密度随时间的变化

Figure  7.  Variation of k_eff and Fuel Uniform Density with Time for Vertical Impact

下载: 全尺寸图片 幻灯片

在变形过程中，对临界影响最大的就是燃料棒的形变。除特征值变化外，图7还展示了垂直跌落时燃料区均匀密度随时间的变化。可以看出，经过碰撞后的挤压，燃料棒的密度由初始的9.9 g/cm³增加为均匀密度下的10.8 g/cm³。图8为燃料区末尾时刻燃料区的非均匀密度分布情况，密度取值范围为9.9~14.0 g/cm³，85根燃料棒互相之间差别不大，但每根燃料棒的轴向的非均匀密度效应显著，与地面接触处的密度远远高于远离碰撞点区域的密度，临界计算按燃料棒均匀密度处理会带来较大的误差，本文对碰撞末尾时刻燃料棒均匀密度变化和非均匀密度变化的k_eff变化情况进行了对比分析，如表3所示，不考虑密度变化的情况是指仅仅考虑变形导致的k_eff的变化，均匀密度变化是在体积变形的条件下，按照质量守恒对燃料等的密度进行统一调整之后再进行临界计算的k_eff变化，非均匀密度变化是在基于每个网格的体积变形的条件下，每一个网格按照质量守恒定律推导出每个网格新的密度，精细地考虑了燃料各个位置的密度情况下的k_eff的计算，可以看出，单纯碰撞引起的k_eff也是增加，叠加了密度效应后引入更大的反应性，并且非均匀密度下反应性的引入比均匀密度下的提升约10%。

图  8  垂直碰撞末尾时刻燃料区非均匀密度分布

Figure  8.  Non-uniform Density Distribution of Fuel Change at the End of Vertical Impact

下载: 全尺寸图片 幻灯片

表  3  垂直碰撞末尾时刻不同处理方式下特征值变化的对比

Table  3.  Comparison of Different k_eff Changes at the End of Vertical Impact

参数	不考虑密度变化	均匀密度变化	非均匀密度变化
Δkeff/pcm	+140	+7993	+8800

下载: 导出CSV 

| 显示表格

3.3   45°倾角跌落分析

相较于空间堆垂直跌落情况，倾斜跌落更加符合实际，ABAQUS模拟最终位移结果如图9所示。可以看出，应力分布存在明显的挤压现象，空间堆外层在与地面碰撞中出现了非常大的变形，图10和图11展示了反应堆跌落前后的结构形状和归一化中子注量率分布的变化，可以看出，碰撞后中子注量率峰值提升了8%。图12也反映出特征值在碰撞前段部分并没有发生太大变化，末尾处燃料发生变形时骤升，MCNP-UM和DAG-OpenMC计算值与文献[12]结果较为一致，与垂直跌落分析类似，差异主要还是来自于文献数据描述的不确定性以及内部UM和表面UM处理方法的不同，k_eff从0.865提升为0.895，增加了约3%。

图  9  45°倾角跌落位移云图

Figure  9.  Displacement Distribution after 45° Impact

下载: 全尺寸图片 幻灯片

图  10  45°倾角跌落碰撞前后反应堆几何结构图

Figure  10.  Reactor Geometry Structure after 45° Impact

下载: 全尺寸图片 幻灯片

图  11  45°倾角跌落碰撞前后中子注量率归一化分布

Figure  11.  Normalized Neutron Flux Rate Distribution before and after 45° Impact

下载: 全尺寸图片 幻灯片

图  12  45°倾角跌落过程中特征值和燃料区均匀密度变化情况

Figure  12.  Variation of k_eff and Fuel Uniform Density with Time for 45° Impact

下载: 全尺寸图片 幻灯片

同时，图12也反映了45°倾角跌落过程中燃料区均匀密度随时间的变化情况。经过碰撞后的挤压，燃料棒的密度由初始的9.9 g/cm³增加为均匀密度下的10.8 g/cm³。图13为燃料区末尾时刻燃料区的非均匀密度分布情况，密度畸变主要集中在与地面接触的部位，未与地面接触的上部燃料密度没有发生变化，非均匀密度情况下燃料密度取值为9.8~14.0 g/cm³，表4所示为碰撞末尾时刻燃料棒均匀密度变化和非均匀密度变化k_eff变化情况的对比，可以看出，非均匀密度下k_eff的增加会提升约20%，与垂直撞击不同，45°倾角跌落引起的是负反馈效应。

图  13  45°倾角跌落末尾时刻燃料区非均匀密度分布

Figure  13.  Non-uniform Fuel Density Distribution at the End of 45° Impact

下载: 全尺寸图片 幻灯片

表  4  45°倾角跌落末尾时刻不同处理方式特征值对比

Table  4.  Comparison of Different k_eff Changes at the End of 45° Impact

参数	不考虑密度	均匀密度变化	非均匀密度变化
Δkeff /pcm	−2276	+3090	+3690

下载: 导出CSV 

| 显示表格

4.   结　论

本文以空间核反应堆系统高速撞击地面的极端变形为例，针对NaK冷却空间堆进行了冷态（控制鼓吸收体面向堆芯）、垂直跌落和45°倾角跌落碰撞的核临界安全分析，基于均匀和非均匀密度方式，建立了MCNP、DAG-OpenMC与ABAQUS耦合分析模型，研究了微型反应堆动力学冲击下的反应性变化，具体结论如下：

（1）在计算未变形情况下，UM蒙特卡罗与传统CSG几何蒙特卡罗的计算精度相当；采用UM蒙特卡罗方法与有限元技术耦合可以进行大变形条件下反应堆的模拟，而采用CSG几何的蒙特卡罗方法则难以模拟。

（2）大变形下空间堆特征值存在一定程度的增加，以NaK冷却空间堆100 m/s速度撞击混凝土为例，垂直撞击引起的k_eff增加约8%，45°撞击引起的k_eff增加约3%，与文献符合较好；内部UM处理能体现材料的非均匀密度效应，燃料非均匀密度变化条件下2种场景的k_eff增加分别提升约10%和20%。

本文的研究对于空间堆等此类微型反应堆的安全设计提出了新的挑战，由于空间堆等微型反应堆碰撞会造成堆芯的较大反应性的引入，需要在设计上对此类事故进行防范，研究成果可以为后续空间堆发射、移动微堆特殊场景的安全评估和新质堆研发及复杂几何辐射场生成提供有力的理论支撑。