核反应堆结构材料网格精细活化分析研究

李雷鸣 江书益 陈珍平 苑旭东 杨 超 谢金森 于 涛

1（南华大学核科学技术学院 衡阳421001）

2（湖南省数字化反应堆工程技术研究中心 衡阳421001）

1 结构材料活化分析方法

在反应堆运行过程中，堆内结构材料受到中子的辐照而发生活化反应，使其产生大量放射性核素，是反应堆屏蔽设计、检修换料方案、退役策略中都需要考虑的重要内容［1］。由于反应堆技术的发展，各式的先进反应堆具有更加复杂的结构以及材料布置与中子能谱，导致中子活化分析也更加复杂，所以有必要开展核反应堆结构材料精细活化分析方法研究［2］。传统的栅元活化方法将每个栅元的中子注量率取平均值进行活化，无法解决由于栅元内部结构材料活化的非均匀空间效应而引起的活化源项及活化γ源空间分布不够精确的问题，进而无法对活化γ辐射场进行精确评价［3］，对比见图1。

图1 核反应堆结构材料活化分析方法对比Fig.1 Comparison of structural material activation analysis methods

为了提高核反应堆结构材料活化分析精度，本文发展了基于网格活化的核反应堆结构材料精细活化分析方法。所开发程序在支持直角网格与柱坐标网格的基础上支持多个非均匀网格拼接，使网格划分更加灵活的同时也提升了网格的上限。本文基于严格两步法（Rigorous 2-step method，R2S），基于蒙特卡罗核粒子输运程序系统（Monte Carlo NParticle Transport Code System，MCNP）［4］与活化程序FISPACT［5］耦合，开发了由“中子输运”到“材料活化”再到“活化剂量”进行耦合的反应堆结构材料精细化活化计算分析程序（MCNP-FISPACT coupled Mesh-based Activation code，MCFisMA）。

2 基于R2S的网格活化分析方法

2.1 严格两步法

本文基于严格两步法建立核反应堆结构材料网格精细活化分析，其主要包括三步：中子输运网格计算、网格活化源项计算和衰变光子剂量计算，流程如图2所示。

图2 基于R2S的网格活化分析方法Fig.2 Mesh-based activation analysis method based on rigorous 2-step method

1）中子输运网格计算：建立核反应堆结构材料蒙特卡罗中子输运计算模型，对结构材料区域构建网格计算模型，确定中子能群结构，基于MCNP 程序开展蒙特卡罗中子输运网格计算，得到结构材料区域的网格中子注量率和网格中子能谱。

2）网格活化源项计算：建立针对结构材料网格模型的网格活化计算模型，根据中子输运计算得到的网格中子注量率和网格中子能谱，基于FISPACT开展网格活化计算，得到网格活化源项精细空间分布，主要包括网格主要活化核素及贡献程度、活化光子源强、活化光子能谱等。

3）活化光子剂量计算：根据2）中得到的网格活化衰变光子源分布，基于MCNP 进行二次开发建立网格衰变光子源抽样方法，构建衰变光子剂量计算蒙特卡罗光子输运计算模型，调用MCNP 进行网格活化源项衰变光子活化剂量场计算。

在网格活化分析方法中，有两大关键问题需解决：1）网格材料等效均匀化计算。针对结构材料网格模型中存在的非均匀材料网格，需要确定其准确材料成分才可实现后续的网格活化计算；2）网格活化光子源精确抽样。对于结构材料的网格活化衰变γ源，需要建立网格衰变光子源抽样方法以实现活化γ剂量场精确计算。

2.2 网格材料等效均匀化计算

核反应堆结构材料模型的几何结构通常较为复杂，几何边界多为异形结构。因此，基于网格模型开展结构材料活化分析时，在建立结构材料区域的网格模型时，网格模型边界可能不一定与结构材料栅元的外边界完美契合，存在一个网格包含多个非规则材料区域情况（图3），这将导致一个网格由多种不同的材料组成，此时则无法直接得到其网格材料核素成分，进而无法直接开展网格活化计算。因此，需要对网格材料进行均匀化处理，以得到每个网格中精确的等效均匀材料，用于后续网格活化计算。具体均匀化如图3所示。

图3 网格材料等效均匀化计算方法Fig.3 Calculation method of equivalent homogenization for mesh material compositions

在计算网格等效均匀化材料时，以当前网格中各种材料的体积比例作为权重因子，通过体积加权，可以得到网格等效均匀化材料中各核素的质量密度或核子密度。网格等效均匀化材料计算公式如下：

式中：mmesh，j与ρmesh，j分别为网格内第j种核素的质量与密度；i为网格内第i种材料；n为当前网格内总材料数；mi，j与ρi，j为当前网格内第i种材料中第j种核素的质量与密度；Vi为当前网格内第i种材料的体积；Vmesh为当前网格的总体积。

网格等效均匀化材料计算流程如图4所示。为了获得每个网格中各核素成分的体积加权因子，本文采用源粒子轨迹统计法实现体积加权因子的计算。针对每个网格，在网格区域内进行一定数目源粒子均匀抽样，并进行蒙特卡罗粒子输运，统计粒子在网格内不同材料区域内的分布情况，进而得到网格内不同材料的体积比例，最后按照式（1）进行网格等效均匀材料的计算。

图4 网格等效均匀材料的计算流程Fig.4 Calculation flow of equivalent homogenization material for mesh

2.3 网格活化γ源粒子抽样

MCFisMA程序进行网格活化计算之后，将产生精细分布的网格活化光子源，通常具有网格源数目较多（102~106个网格源区域）、几何分布复杂（直角网格或扇形网格）等特点，此时MCNP 程序自带的源定义功能将无法满足网格源定义要求。因此，本文对MCNP 程序源粒子抽样模块进行二次开发，以实现以大规模复杂网格γ 源抽样，进而实现网格活化γ剂量场的精确计算。

网格γ源粒子抽样时，主要考虑网格抽样概率、空间位置抽样、运动方向抽样和能量抽样。其抽样流程如下：首先，根据网格的衰变光子源强与整个模型总光子源强比例，得到每个网格的归一化抽样概率，根据抽样概率可抽样并确定源粒子所在的网格；其次，根据图5的流程进行源粒子坐标、运动方向以及能量的抽样，进而实现网格光子源抽样。

3 数值验证

3.1 ITER停堆剂量率基准验证

3.1.1 ITER基准模型介绍

为了验证本文方法及MCFisMA 程序的正确性，将其应用到ITER 国际组织发布停堆剂量（Shutdown Dose Rate，SDR）基准题中。ITER 例题被广泛应用于验证活化程序的正确性，模型结构如图6 所示。在基准题中，由于模型尺度较大且屏蔽材料以铁水混合物为主，存在较为严重的深穿透效应。在蒙特卡罗输运计算过程中采用了重要性减方差方法，首先将屏蔽体分为大小相等的圆环，再将重要性从模型中子源端开始依次增加，首先尝试以指数为倍数增加，再根据输出的tracks entering 和population信息进行调整，使得进入每一个细分栅元的粒子数基本相同，即可以开始中子输运计算。

图6 ITER停堆剂量率基准题Fig.6 ITER shutdown dose rate benchmark

3.1.2 停堆剂量率结果分析

剂量率的计数区域为图6中的4个计数栅元，在进行活化计算时，首先将活化区域划分为圆柱网格，分别在r、θ、z坐标轴上划分10、10、3 个网格，总计300 个圆柱网格进行活化计算。使用14 MeV 的聚变中子源对基准题模型进行辐照，采用SA2-ITER的辐照方案。计算辐照后冷却106s时4个计数栅元的停堆剂量率，计算结果如图7所示，其中核能与应用 实 验 室（Nuclear Energy and Application Laboratory，NEAL）为本文结果。

图7 不同国际机构ITER停堆剂量率对比Fig.7 Comparison of ITER SDR benchmark of different international organizations

图7给出了由本文方法得到的停堆剂量率计算结果与国际上其他机构计算值［3］对比情况。计算结果显示，本文结果趋势吻合较好，数值上与英国卡拉姆聚变能研究中心（Culham Centre for Fusion Energy，CCFE）和日本原子能机构（Japan Atomic Energy Agency，JAEA）最为接近，其中CCFE 为ITER 总结报告的参考单位之一。其细微区别来自于截面库以及活化计算方法的区别，以及不同网格划分方法也会对结果产生一定偏差。基于以上结论，可以证明本文方法和程序的正确性。

3.2 NUREG/CR-6115压水堆精细活化分析

3.2.1 NUREG/CR-6115压水堆模型介绍

相比于传统栅元活化程序，MCFisMA程序可以得到更为精细的结构材料活化光子源分布。为了证明本文方法及程序的精细程度和适用性，将其应用到NUREG/CR-6115压水堆中，其模型见图8。

图8 NUREG/CR-6115压水堆模型Fig.8 The NUREG/CR-6115 PWR reactor model

对NUREG/CR-6115 压水堆主要结构部件如吊篮、围板、热屏蔽层、压力容器以及压力容器内衬进行了网格活化精细计算。其中，网格活化模型使用柱坐标网格，角向划分30个网格，轴向划分128个网格。如图9 所示，MCFisMA 程序给出了堆芯围板、堆芯吊篮、热屏蔽层、压力容器等结构材料区域的活化光子源的空间高分辨率精细分布。数据结果表明，活化光子源的源强随其距堆芯距离的增加呈递减趋势，部件的活化程度也随距堆芯的距离而降低。而压力容器内衬的体积较小，故其衰变光子源的源强远低于其他结构，四个部件的最大光子源强对比见图10，其中横坐标中的1、2、3、4 依次代表堆芯围板、堆芯吊篮、热屏蔽层、压力容器。针对同一结构材料而言，在活性区高度的中心区域，由于该区域中子注量率较高，材料活化程度也较高，衰变光子源的源强也较高。

图9 NUREG/CR-6115压水堆结构部件活化光子源项精细分布 (a)堆芯围板，(b)堆芯吊篮，(c)热屏蔽层，(d)压力容器Fig.9 High-resolution distribution of activation photon sources for NUREG/CR-6115 structural components(a)Baffle,(b)Core barrel,(c)Thermal shield,(d)Pressure vessel

图10 NUREG/CR-6115压水堆结构部件最大光子源强对比Fig.10 Maximum photon source strength comparison of NUREG/CR-6115 PWR structural components

3.2.2 结构材料活化剂量场分布计算

基于MCFisMA 程序分析了不同网格精细程度对NUREG/CR-6115 压水堆结构材料活化剂量分析的影响。主要对堆芯围板、吊篮、热屏和压力容器等结构材料进行不同精细程度的轴向网格划分，分析了不同轴向网格数目下活化剂量场分布变化情况（图11）。随着轴向网格数目（图11 中z表示网格数目）的变化，剂量场分布也呈现变化明显。轴向网格的数目较少时，整个剂量场分布相对均匀，同时也存在较大的误差。随着网格数目增加，剂量场分布曲线逐渐收敛，逼近精确值。其主要原因在于轴向方向的结构布置不均匀，尽管在堆芯上下方布置有反射层，但中子注量率分布仍不均匀。在轴向上，靠近堆芯活性区高度区域中子注量率较大，因此随着网格数目的增加，靠近活性区高度区域内结构材料活化源项结果更加精确，使得活化γ 光子所致剂量分布将逐渐收敛于真实值。

图11 不同轴向网格数目下的剂量场分布Fig.11 Dose distribution under different axial mesh numbers

同时，为了量化分析网格活化方法相比传统栅元活化方法对计算结果的影响，选取了NUREG/CR-6115模型轴向方向的最低点、中间点和最高点进行了不同网格数目下剂量场对比分析，如表1 所示。在网格数为1 的情况下，可将其视为传统栅元活化计算结果。与栅元活化结果相比，随着网格数目增加，剂量相对误差逐渐降低，在网格数达到8之后相对误差将显著降低，当轴向划分为16 个网格时（网格边长24 cm），得到较为理想的剂量率分布结果（综合考虑结果精度与计算时长，因为随着网格数提升计算速度降低）。

表1 不同网格数目(N)下轴向位置剂量率的相对误差(%)Table 1 Relative error of dose rate under different mesh numbers(%)

4 结语

本文针对核反应堆结构材料活化的问题，开展了基于两步法的网格模型精细活化分析方法研究，考虑了中子注量率在空间中的非均匀分布，实现对结构材料活化源项和活化γ源的精细网格空间分布计算，进而实现活化源项精细分析和活化γ 辐射场的精确评价。同时，基于MCNP与FISPACT耦合开发了全自动耦合网格活化分析程序MCFisMA，通过ITER停堆剂量率基准题验证程序的正确性，并基于NUREG/CR-6115 压水堆模型实现了结构材料网格精细活化分析。通过网格活化分析数值结果可知，不同精细程度网格模型对活化光子主要贡献核素、活化光子能谱、活化光子源强均会产生较大影响，尤其对活化剂量场分布影响显著。相比于传统栅元活化分析方法，本文网格活化分析方法可以获得活化核素、光子源等源项的精细空间分布，可评估核反应堆复杂结构部件中任意空间位置处的源项分布情况，为反应堆屏蔽设计、检修换料方案、退役策略提供理论与数据支撑。

作者贡献声明李雷鸣：主要负责网格活化计算理论研究、程序开发、起草文章等工作；江书益：主要负责程序数值基准验证研究；陈珍平：负责程序在压水堆中应用及网格活化分析研究，提供研究经费，对文章进行审阅；苑旭东：负责网格活化源方法及源抽样程序开发；杨超：负责对网格源程序验证研究；谢金森：负责蒙卡粒子输运计算理论理论指导；于涛：负责蒙卡粒子输运减方差理论指导、提供研究经费及超算平台，对文章内容进行审阅。