核反应堆屏蔽计算堆芯中子源强生成方法研究

张平逊 张 斌 陈义学

（华北电力大学 核科学与工程学院 北京 102206）

离散纵标法（SN）是粒子输运方程求解中重要的确定论方法之一，在解决深穿透屏蔽问题中具有显著优势［1］。源强作为重要的输入参数，计算精度直接影响屏蔽计算的可靠性。源强受几何模型、燃耗与功率分布的影响［2］，计算结果不确定性［3-4］较大。全堆芯Pin-by-Pin源强处理可较准确地描述堆芯源强分布，但几何建模复杂繁琐，实际应用并不理想。针对堆芯中子源强计算方法，国内外开展了大量研究。Haghighat［5］研究发现，当源强计算考虑堆芯不同裂变核素与采用不同核素处理方法时，将对结果最高引入10%以上的偏差。西屋公司综合考虑源强计算相关影响因素，为DORT［6］等离散纵标程序开发了SORCERY前端源项处理程序，但受计算机内存限制，无法处理大规模三维问题［7］。上海核工程研究设计院完善了TORT程序的源项生成功能，提高程序的应用范围［8］。总体上，国内外源强计算方法研究已较为成熟，为获取精确的堆芯源强分布，针对大规模屏蔽问题采用Pin-by-Pin源强计算时，需对燃料栅元进行大量冗余几何建模，计算效率存在较大挑战。

本文研究基于多维离散纵标输运计算程序ARES［9］进行展开，分析源强计算方法对屏蔽计算精度的影响。以平均源强计算方法为基础，构造多权重源强网格映射算法，避免复杂Pin-by-Pin几何建模，提高源强计算效率与精度，降低程序内存限制。

1 理论模型

1.1 中子源强计算方法

源强计算作为SN方法屏蔽计算的前处理过程，输运计算前需进行几何网格的源强赋值，基本公式为：

式中：χig为固定源能谱；Ci为功率源强转换因子；Pi为燃料组件功率；i为堆芯组件编号；g为能群编号。针对单一裂变核素进行源强计算，功率源强转换因子定义为：

式中：υ为核素每次裂变平均释放的中子数；Er为每次裂变平均释放的能量。由于组件燃耗变化，如裂变核素钚的产生与沉积将对堆芯源强产生重要影响［10-11］。考虑各裂变核素对源强的贡献，采用截面平均方法对固定源能谱与功率源强转换因子进行加权平均［5］：

1.2 多权重源强网格映射算法

1.2.1 中子价值理论

根据中子守恒原理，屏蔽问题中不含裂变源的稳态固定源中子输运方程表示为［12］：

式中：从左往右依次为泄漏项、移出项、散射源项和外中子源项为中子角通量密度为宏观总截面为宏观散射截面分别为空间、能量和角度变量。

式中：H*为H的共轭算子。根据价值守恒理论推导共轭中子输运方程，详细过程见参考文献［1］：

对中子输运方程与共轭输运方程进行简写：

分别将式（8）和式（9）与ψ*和ψ作内积，并将两式相减，边界条件为真空边界，则有：

若函数Σd(E)为中子与某探测器的响应值，则共轭函数表示探测器对此中子的响应值，即中子价值。屏蔽计算中，通过一次共轭计算即可获得堆内各燃料组件对压力容器处探测器响应的贡献。

1.2.2 多权重源强网格映射算法

反应堆各组件源强对压力容器处中子注量率的贡献不同［13］，采用共轭计算可获得堆芯共轭通量密度分布，共轭通量密度大的区域其中子价值高，产生的中子对压力容器处中子注量率的影响较大，几何模型需以燃料栅元为单元进行Pin-by-Pin精细建模。

在压力容器处构建共轭源，经共轭输运可得堆芯源强区域中子价值分布。如图1所示，燃料组件由外向内中子价值迅速降低，外围组件产生的中子对压力容器处中子注量率起主要贡献。多权重源强网格映射算法充分考虑堆芯组件中子价值分布，对中子价值大的外围组件进行精细源强计算。源强网格与几何网格在相同笛卡尔直角坐标系下处理，与几何网格划分不同，源强网格为程序自动处理计算的区域坐标。几何网格源强赋值流程如图2所示，主要分为以下三步：

图2 几何网格源强赋值流程Fig.2 Flow chart of geometric grid source assignment

1）由中子价值对功率分布进行近似：外围两层组件中子价值分布为0.1~1.0，对输运计算产生显著影响，组件使用Pin-by-Pin计算，内层组件的中子价值小于0.1，对输运结果影响较小，可采用均匀近似，经用户定义的几何输入卡获得不同层级的源强区域边界与功率分布；

2）根据各组件燃耗信息，使用截面平均方法对源强裂变谱与功率源强转换因子进行加权平均，结合源区边界与功率获得区域源强；

3）遍历源区各几何网格，采用映射算法实现区域源强向几何网格的源强赋值，其中非源区网格源强直接赋值为零，节约时间成本。

源强区域划分为计算机基于功率分布自动处理的虚拟区域，实际模型并未进行划分，且源区边界与网格边界位置不同，无法直接进行源强赋值。为实现源强赋值，保证总源强守恒，采用体积权重法对几何网格进行分割，并对划分的子网格分别进行源强计算。其中有：

式中：Smesh为几何网格源强；Vmesh为几何网格体积；Si为按源强区域边界划分的子网格源强；Vi为子网格体积；M为划分的子网格数目。如图3所示，堆芯某几何网格跨多个源强区域，按源区边界将网格进行分割。分别对分割产生的8个子网格区域进行源强计算，采用式（13）进行几何网格的源强赋值求和。依次对各燃料组件进行网格源强映射计算，实现堆芯几何网格源强赋值。

图3 几何网格源强映射方法 (a) 几何网格分割，(b) 子网格源强计算，(c) 几何网格源强赋值Fig.3 Geometric grid source mapping method (a) Geometric grid segmentation, (b) Sub-grid source calculation, (c) Geometric grid source assignment

2 基准验证

2.1 基准题描述

NUREG/CR-6115基准题［14］中标准堆芯燃料装载方案为典型压水堆模型，包含204个燃料组件，总功率为2 527.73 MW，该基准题可用于屏蔽程序的评估与分析。基准题模型堆芯外围组件径向功率分布梯度较大，且轴向功率峰因子对源强结果影响明显，使用多权重源强网格映射算法与平均源强算法分别进行计算，分析多权重源强网格映射算法的计算精度与效率。图4为NUREG/CR-6115基准题的几何模型，反应堆由堆芯向外包含吊篮、下降段、热屏、水隙、压力容器、空腔和混凝土，压力容器内半径为219.075 cm，厚21.59 cm。反应堆内其他结构尺寸和材料成分在NUREG/CR-6115基准报告中详细给出。

图4 NUREG/CR-6115基准题几何模型1-围板，2-吊篮，3-内部水隙，4-热屏，5-压力容器内衬，6-空隙，7-压力容器隔热层，8-混凝土内衬Fig.4 Geometric diagram of the NUREG/CR-6115 benchmark model 1-Core baffle, 2-Core barrel, 3-Inner inlet water gap,4-Thermal shield, 5-Reactor pressure vessel liner, 6-Air gap, 7-Pressure vessel insulation, 8-Biological shield liner

基于KASHIL-E70多群截面数据库［15］进行截面处理，采用离散纵标屏蔽程序ARES进行三维输运计算。网格划分为284×284×161，空间离散格式选择指数定向θ权重差分方法（Exponential Directional Weighted difference method，EDW），各向异性散射截面使用P3阶勒让德多项式进行展开近似，求积组采用勒让德-切比雪夫求积组（PNTN），求积阶数设置为S8，迭代收敛精度为1×10-3。

2.2 结果分析

堆芯源区以燃料组件为基本单元进行几何建模，使用多权重源强网格映射算法与组件平均源强计算方法分别进行源强计算，输运参数设置相同。图5为两种源强计算方法下三维堆芯源强分布，多权重源强网格算法轴向源强分布由中间位置向两端递减，径向外围组件源强呈梯形分布，相比于平均源强计算方法能更准确描述堆芯源强。选取压力容器内表面轴向峰值处与焊缝处探测点位置的中子能谱分布进行分析，能谱分布结果如图6所示。压力容器轴向峰值处中子能谱分布趋势总体一致，吻合程度相对较好，但焊缝处中子能谱分布吻合相对较差，平均源强计算方法输运结果偏大。偏差产生的主要原因为轴向功率峰因子对不同轴向区域快中子注量率的影响不同，造成输出结果稳定性较差，图7为外围组件轴向功率分布。

图5 NUREG/CR-6115模型源强分布 (a) 平均源强计算方法，(b) 多权重源强网格映射算法Fig.5 Source distribution of the NUREG/CR-6115 model(a) Average source calculation algorithm, (b) Multi-weight source mesh mapping algorithm

图6 压力容器内表面中子能谱分布 (a) 压力容器轴向峰值处，(b) 压力容器焊缝处Fig.6 Neutron spectrum distribution on internal surface of pressure vessel (a) Axial peak location, (b) RPV lower weld location

图7 NUREG/CR-6115模型轴向归一化相对功率分布Fig.7 Axial-normalized relative power distribution of the NUREG/CR-6115 model

为比较两种源强计算方法的优劣，对总体输运计算结果进行分析。选取压力容器不同轴向位置随方位角变化（0°~45°）的探测点统计快中子注量率（E＞1.0 MeV）输运结果，按式（14）计算与基准题报告中MCNP和DORT参考值的相对误差均方根：

式中：φcalc，i为输运计算值；φref，i为基准参考值；N为选取的探测点数目。

表1为NUREG/CR-6115模型压力容器处不同轴向位置快中子注量率ERMS结果，总体上多权重源强网格映射算法与平均源强计算方法相比ERMS偏小。其中多权重源强网格映射算法输运结果与DORT参考值相比ERMS最大约为14.4%，与MCNP参考值相比最大约为5.7%；平均源强计算方法输运结果与DORT参考值相比ERMS最大约为29.7%，与MCNP参考值相比最大约为21.2%。由于平均源强计算方法将各组件源强进行平均，造成压力容器不同轴向高度位置快中子注量率ERMS不同，焊缝处ERMS比轴向峰值处增大近一倍。相比于平均源强计算方法，多权重源强网格映射算法输运结果更精确，其中压力容器焊缝处计算结果与MCNP参考值相比ERMS降低约18.46%。

表1 反应堆压力容器快中子注量率（E＞1.0 MeV）相对误差均方根Table 1 Root mean square (RMS) of the relative error for fast neutron fluence in reactor pressure vessel (E＞1.0 MeV)

多权重源强映射方法建模方便，使用灵活的网格映射技术，可实现任意位置源强赋值。以外围三层组件Pin-by-Pin输运结果为参考值，选择压力容器处探测点位置，分别对比不同分层下快中子注量率分布结果。表2为源强映射方法下不同层级数源强相对误差分布，能量E＞1.0 MeV时，平均源强计算方法与外围三层组件Pin-by-Pin快中子注量率相对偏差的平均值为20.68%，能量E＞0.1 MeV时相对偏差为19.7%。外围一层组件与外围两层组件平均相对功率偏差小于0.5%，其中外围两层组件输运计算结果与参考值基本吻合，满足计算精度要求。

表2 中子源强分布影响Table 2 Influence of neutron source distribution

3 结语

本文针对屏蔽计算中复杂源强处理问题，提出了多权重源强网格映射算法。通过对几何权重较大的组件进行精细计算，并采用体积权重法将区域源强向几何网格进行映射，保证总源强守恒。避免复杂几何建模，提高计算效率，并有效降低计算机内存限制。数值计算表明，多权重源强网格映射算法相比于平均源强计算方法，提高了输运计算的精度，且保障轴向位置计算结果的稳定性。另外，光子输运在屏蔽计算中占有重要作用，光子源强将是未来研究的工作重点。

作者贡献声明张平逊负责文献调研、程序功能开发与测试，以及文章撰写和修订；张斌负责相关文献资料的搜集与整理并对文章进行修订；陈义学提供程序研发平台。