节块法与蒙特卡罗方法耦合计算研究

王新哲，喻 宏，徐 李，胡 赟，杨晓燕

（中国原子能科学研究院快堆研究设计所，北京 102413）

节块法与蒙特卡罗方法耦合计算研究

王新哲，喻 宏，徐 李，胡 赟，杨晓燕

（中国原子能科学研究院快堆研究设计所，北京 102413）

普通节块法无法在计算中获得不同组件内精细中子通量密度分布的信息。本文提出一种利用入射角通量将节块法与蒙特卡罗方法相耦合的方法（节块－蒙卡入射角通量耦合方法），并编制了计算程序进行验证。结果表明：本文计算结果与参考值相符，节块－蒙卡入射角通量耦合方法适用于局部特定位置精细中子通量密度等参数的计算，计算效率高，计算结果准确。

节块法；蒙特卡罗方法；耦合；入射角通量；精细中子通量密度

堆芯精细中子通量密度和功率分布的计算是反应堆中子学分析软件发展的方向之一，准确获得堆芯精细中子通量密度和功率分布，可有效减少设计中的保守量，提高反应堆的经济性。普通节块法在计算中无法获得不同组件内精细中子通量密度分布等参数的信息。现阶段，精细中子通量密度与功率分布计算的主要方法有功率重构方法［1］、蒙特卡罗方法［2］及特征线方法［3］等。这些方法在现有条件下皆存在一定局限性，如功率重构方法对于复杂结构组件计算精度差［4］，蒙特卡罗方法与特征线方法的计算速度较慢。

本文通过分析节块法与蒙特卡罗方法的特点，提出一种利用入射角通量将二者相耦合的方法——节块－蒙卡入射角通量耦合方法，对任意几何组件的精细中子通量密度与功率分布进行计算。

1 理论模型

1.1 节块－蒙卡入射角通量耦合方法描述

节块－蒙卡入射角通量耦合方法的基本思想是，将求解中子输运方程的两种思路——确定论方法与蒙特卡罗方法通过边界流数据传递的方式耦合在一起，从而发挥二者各自的优势。

稳态中子输运方程为：

其中：Φ为中子通量密度；r为中子空间位置向量；Ω、Ω′分别为碰撞前后中子运动方向的单位向量；E、E′分别为碰撞前后中子的能量；Σt为材料总宏观截面；Σs为材料散射截面；f为散射函数；Qf为裂变源项；S为外源项。

确定论方法的基本思想是离散，中子输运方程中的未知参数中子通量密度Φ（r，E，Ω）是空间位置、能量与飞行角度的函数，这些自变量在物理上都是连续的。为了将连续的中子输运方程转换为有限个可解的数学物理方程，必须对这些连续参数进行离散。现有确定论方法对于几何空间的离散，大多采用规则网格划分，不能进行精细几何描述，对于结构复杂的组件设计有一定计算误差。

蒙特卡罗方法能进行精细几何描述，问题适应能力强，但由于其需对大量粒子进行统计，因此全堆计算需较长时间。

在某些计算中只关注堆内特定位置处的精细分布，上述两种方法很难快速、准确地给出计算结果。此时，可利用将二者耦合的思路，具体计算流程如下：1）利用节块法程序进行快速全堆计算，得出中子通量密度分布和有效增殖因数等信息；2）将获得的信息转换成所求关注位置处组件的入射角通量；3）利用获得的入射角通量信息产生蒙特卡罗程序的面源文件，同时按照精细几何进行组件描述，进行蒙特卡罗面源计算。

通过上述3步，即可以较低的计算代价获得特定组件的精细几何描述计算参数。

1.2 入射角通量分布计算

本文采用六角形节块程序NAS［5］作为节块法程序。NAS由中国原子能科学研究院快堆研究设计所开发，能按照扩散与输运两种方法进行三维堆芯计算，具有微扰计算、时空动力学计算、燃耗计算、堆芯燃耗管理及优化等功能。NAS在径向上按照组件划分节块，轴向则可根据堆芯布置适当进行划分。本文堆芯计算中，采用NAS输运模块进行，表面通量采用双球谐DP3近似展开，节块的求解采用了横向积分方法，它将三维输运方程转化为4个一维输运方程，并采用一维离散纵标（SN）方法进行求解。

NAS通过堆芯计算，可给出节块各面上的入射通量Φin、入射流Jin、二阶入射通量Φin2及三阶入射通量Φin3，节块间的耦合关系为：

其中：j′和j分别为临近的两个节块；ψij为j节块i方向的入射角通量；μi和ωi分别为SN方法中高斯－勒让德求积组的离散方向余弦与求积权重，其中μi转换到［0，1］区间，∑iωi＝1；Pn为n阶勒让德多项式。在现版本NAS中，角度的离散只考虑了与节块表面法线夹角一个方向。

据此可得到入射角通量的表达式：

其中：

通过式（3）与式（4）即可求得各面上分能群入射中子角通量ψ（Si，Ej，μl），其中，Si为第i面，Ej为第j能群，μl为第l离散角度。

要获得蒙特卡罗程序所需的面源分布，需将ψ（Si，Ej，μl）转换为面入射粒子数，按照式（5）进行转换。

其中：N（Si，Ej，μl）为面入射粒子数；Ai为入射面面积。

对于总计算粒子数NT，各面上、各能群、各角度的入射中子数N′（Si，Ej，μl）有：

1.3 蒙特卡罗面源抽样方法

本文中蒙特卡罗程序采用MCNP［6］，MCNP可读取面源文件RSSA进行输运计算。本文根据式（6）中的中子数分布通过蒙特卡罗方法进行抽样，自动生成该文件。面源文件需提供每个中子的编号、时间、权重、空间坐标、飞行方向、与入射面夹角和能量等信息，其中需要抽样的信息为空间坐标、飞行方向和能量。

由于NAS尚不能提供入射中子面分布信息，故在抽样过程中认为每个入射面上中子的分布为均匀分布。为方便抽样，抽样方法为先在一个标准矩形面上进行抽样，然后再将粒子坐标转换到相应的面上。

对于角度分布，认为在每个离散角度内，中子飞行方向均匀分布，即：

其中，f（φ（μ））为飞行方向为μ的中子通量密度的概率密度。

则入射粒子与面夹角的概率密度函数为：

其中，C为归一化常数。

由式（7）、（8）可见，入射粒子关于面夹角的分布为一次分布。抽样过程中也需进行相对面位置的转换，在图1所示编号为1的面上，按照各向同性空间粒子抽样［7］，然后按式（8）以一定概率对其取舍，获得夹角参数sinθ、cosθ、sinφ、cosφ，其中θ为极角，φ为方位角，则由式（9）可得到粒子的飞行方向。

图1 六角形节块面编号示意图Fig.1 Scheme of surface number for hexagonal nodal

其中，U、V、W分别为中子与x、y、z轴夹角的余弦。

按照式（10），对上述夹角进行旋转编号，即可得到粒子飞行方向。

式中，N为面编号。

六角形节块上顶面、下底面的入射中子分布可很容易地抽样得到。

各能群内能量的分布由程序输入给出，可为裂变谱、1／E谱、麦克斯韦谱或均匀分布谱。其中，裂变谱采用乘减抽样方法，1／E谱和麦克斯韦谱采用乘抽样方法。本文中NAS采用171群的NVITAMIN－C库，为方便快堆计算，将最后11群并为一群，每群内能量分布按照NJOY程序手册中用于VITAMIN－E库制作的权重函数给出，并将其中高能部分的聚变谱改为1／E谱［8］。

2 程序开发

按照上述方法，基于NAS与MCNP编制了用于六角形组件计算的MCNTC程序。MCNTC程序流程图如图2所示。

图2 MCNTC程序流程图Fig.2 Flow chart of MCNTC

由图2可见，只需给出NAS输入文件、MCNP组件描述文件及部分控制参数，即可完成程序功能。通过输入可控制抽样的能群划分及群内中子能谱、中子飞行方向离散角度划分等参数。

3 程序验证

采用3个算例对节块－蒙卡入射角通量耦合方法及MCNTC进行验证。

3.1 中国实验快堆首炉堆芯燃料组件内功率分布验证

中国实验快堆（CEFR）是我国第一座钠冷快中子反应堆，设计热功率为65MW，于2010年7月首次达到临界。在初装料时，堆芯装有79盒燃料组件、8个控制棒组件和1个中子源组件，堆芯外面为钢反射层，共有336个反射层组件。反射层外为屏蔽层，装载有230个硼屏蔽层组件［9］。燃料组件每盒有61根燃料棒，燃料芯块为富集度64.4%的UO2。堆芯燃料组件编号如图3所示。

图3 堆芯燃料组件编号示意图Fig.3 Scheme of fuel assembly number for reactor core

按照CEFR物理启动实验测得的临界棒位，MCNTC中采用NAS输运功能进行堆芯计算，keff的计算值为1.004 03，与实验值符合较好。耦合计算中MCNP计算采用燃料组件活性区精细描述，并进行裂变功率统计，每次计算抽样107个面源粒子，棒功率统计不确定度小于0.25%。同时，采用MCNP全堆计算作为参考结果，MCNP全堆计算中，燃料活性区为精细几何描述，燃料其他区域及其他组件采用均匀几何描述，并统计所有燃料棒的裂变功率，有效统计中子数为108，keff的计算值为0.999 91±0.000 07，棒功率统计不确定度小于0.20%。

MCNTC与MCNP皆在Intel Xeon E7－4850上运行。MCNTC采用串行计算，每次程序运行时间约为30min，其中约10min为堆芯临界计算，其余约20min为面源生成及组件精细功率分布计算，79盒燃料组件全部计算约需40CPU小时。MCNP采用10核20线程并行计算，程序运行时间约为2 500CPU小时。因此，MCNTC在全堆功率分布计算上的速度较传统MCNP快数十倍。若只计算个别组件，计算效率还能提升。

MCNTC给出79盒燃料组件中各棒的功率分布，选取1／3堆芯，MCNTC与MCNP的计算结果对比如图4所示。

图4 MCNTC与MCNP计算的组件功率分布结果对比Fig.4 Comparison of power distributions for MCNTC and MCNP calculation

图4中：MAX为棒相对功率最大偏差，按式（11）计算；AVE为棒功率代数平均偏差，按式（12）计算；RMS为棒功率几何平均偏差，按式（13）计算；PDIF为最高功率棒偏差，按式（14）计算。

对计算结果进行统计，MCNTC与MCNP计算结果的代数平均偏差小于0.8%，棒相对功率最大偏差小于1.9%，对79盒组件的统计数据列于表1。

表1 MCNTC与MCNP计算结果比较Table 1 Comparison of calculation results for MCNTC and MCNP

图5 组件内的相对功率分布与相对偏差Fig.5 Relative power distribution and relative deviation of assembly

图5示出采用MCNTC计算得到的2－1、6－29位置组件内相对功率分布及与MCNP计算结果的棒功率相对偏差。由图5可看出，组件功率峰通常出现在燃料的边角位置，且靠近堆芯中央。同时，组件位置越远离堆芯，组件功率峰因子越高。

表2列出部分组件位置MCNTC与NASPIN［10］计算结果的对比。表2中给出的结果均为与参考解的比值，参考解由MCNP计算。由表2可知：对于AVE，NAS－PIN的计算结果较MCNTC的稍好；对于PDIF，MCNTC的计算结果普遍较NAS－PIN的好。

3.2 CEFR核反应率分布验证

CEFR在物理启动实验过程中，进行了核反应率分布测量［11］。实验过程中发现，197Au的俘获相对反应率通过MCNP程序F4卡获得的理论计算值与实验值在反射层区存在较大偏差。经过理论分析认为，造成该偏差的主要原因是197Au的辐射俘获截面在5eV附近存在一很强的共振峰，处于该共振峰峰值处的中子自由程与探测片厚度相近，故必须考虑该空间自屏影响，而之前的理论计算中未考虑该影响，导致在反射层区域计算反应率约是实验值的2倍。

表2 MCNTC与NAS－PIN计算结果对比Table 2 Comparison of calculation results for MCNTC and NAS－PIN

采用MCNTC进行核反应率分布计算。其中，燃料区辐照组件横截面如图6所示。在辐照组件辐照管中部放置一直径为5mm、厚度为30μm的金箔。堆芯区辐照组件位置为2－2、3－3、4－3、5－3、6－4，反射层区域辐照组件位置为7－5、9－6、11－8。燃料区辐照组件抽样108次，反射层区辐照组件抽样2×108次。

不同测量位置下，探测片相对反应率分布如图7所示，实验值与计算值分别以2－2位置归一。由图7可见：计算值与实验值符合较好，最大相对偏差小于5%；197Au的俘获反应率在燃料区基本保持不变，而在反射层区明显上升，这主要是因为该区域有较多中子慢化到共振区能量，共振中子对反应率贡献很大。

图6 燃料区辐照组件横截面Fig.6 Layout of test assembly in fuel zone

图7 不同位置的俘获反应率分布Fig.7 Distribution of capture reaction rate in different positions

3.3 CEFR金箔探测片共振自屏验证实验设计及结果

为进一步验证金箔的共振自屏效应，在2014年CEFR提升功率过程中，设计了探测片测量实验。该实验将7片直径为9mm、不同厚度的金箔叠放在一起，置于辐照组件内，辐照组件位于堆芯第2排反射层8－5位置。于2014年3月14日在CEFR上以0.1%额定功率辐照2h。7片金箔的厚度列于表3。

表3 金箔的厚度Table 3 Thickness of gold foil

辐照后将辐照组件拆解，利用高纯锗谱仪逐个测量探测片中198Au位于411.8keV处的特征峰活度，用于确定197Au的（n，γ）反应率。实验测量的探测片相对比活度及MCNTC计算结果列于表4。表4中的结果均按照活度最小的4号探测片归一。

表4 探测片活度实验值与计算值结果比较Table 4 Comparison of gold foil activities for calculation and test

由表4可知，MCNTC计算结果与实验值符合较好，相对偏差在10%以内。同时，该实验验证了金箔在快堆反射层区域存在显著的共振自屏效应，探测片表层反应率是内部反应率的2倍左右。

4 结论

1）与功率重构方法相比，对于规则组件，本文方法所得结果与其相近，计算速度稍慢。但本文方法对组件几何与材料的适应能力更好，能方便计算各类复杂几何组件。

2）与MCNP相比，在全堆精细功率计算上本方法可得到与之相近的结果，但计算速度提升明显，且组件离堆芯越远，相对效率提升越明显。

3）本文方法可用于某些特殊的应用，如反射层中探测片反应率计算。

［1］ WAGNER M R.Three－dimensional nodal and transport theory methods for hexagonal－Zgeometry［J］.Nuclear Science and Engineering，1989，103（4）：377－391.

［2］ 谢仲生，邓力.中子输运理论数值计算方法［M］.西安：西北工业大学出版社，2005.

［3］ ASKEW J.A characteristics formulation of the neutron transport equation in complicated geometries，AEEW－M 1108［R］.UK：Atomic Energy Authorty，1972.

［4］ FORGET B.A three dimensional heterogeneous coarse mesh transport method for reactor calculations［D］.USA：School of Mechanical Engineering，Georgia Institute of Technology，2006.

［5］ 徐李，马大园，施工，等.快堆三维六角形节块法输运计算研究［J］.原子能科学技术，2013，47（2）：161－165.

XU Li，MA Dayuan，SHI Gong，et al.Research of 3－D hexagonal nodal transport method for fast reactor［J］.Atomic Energy Science and Technology，2013，47（2）：161－165（in Chinese）.

［6］ Team X－M C.MCNP：A general Monte Carlo N－particle transport code，Version 5，LA－UR－03－1987［R］.USA：Los Alamos National Laboratory，2003.

［7］ 许淑艳.蒙特卡罗方法在实验核物理中的应用［M］.北京：原子能出版社，2006.

［8］ MACFARLAN R E.The NJOY nuclear data processing system［R］.USA：Los Alamos National Laboratory，1994.

［9］ 田和春.中国实验快堆堆芯设计说明［R］.北京：中国原子能科学研究院，2002.

［10］曹攀.快堆组件内精细功率分布的计算研究［D］.北京：中国原子能科学研究院，2012.

［11］范振东，陈晓亮，胡定胜，等.中国实验快堆核反应率分布测量研究［J］.原子能科学技术，2013，47（增刊）：107－110.

FAN Zhendong，CHEN Xiaoliang，HU Dingsheng，et al.Experimental research of nuclear reaction rate in China Experimental Fast Reactor［J］.Atomic Energy Science and Technology，2013，47（Suppl.）：107－110（in Chinese）.

［12］曹攀，喻宏，徐李，等.快堆燃料组件内精细功率分布的计算［J］.原子能科学技术，2013，47（增刊）：287－290.

CAO Pan，YU Hong，XU Li，et al.Research on pin power distribution in fuel subassembly of fast reactor［J］.Atomic Energy Science and Technology，2013，47（Suppl.）：287－290（in Chinese）.

Study on Nodal and Monte Carlo Coupling Method

WANG Xin－zhe，YU Hong，XU Li，HU Yun，YANG Xiao－yan
（Department of Fast Reactor Research and Design，China Institute of Atomic Energy，Beijing102413，China）

The conventional nodal method could not calculate the distribution of fine neutron flux density in assemblies.A coupling method，called incident angular flux nodal and Monte Carlo coupling（IFNMC）method，was put forward.The computer code was programmed to verify this method.The calculation results are consistent quite well with the reference data.The IFNMC method is both effective and accurate for fine neutron flux density and other parameter calculations in special assembly.

nodal method；Monte Carlo method；coupling；incident angular flux；fine neutron flux density

TL329.2

A

1000－6931（2015）02－0200－07

10.7538／yzk.2015.49.02.0200

2014－09－16；

2014－12－02

863计划资助项目（2011AA050301）

王新哲（1989—），男，山东东营人，研究实习员，反应堆物理专业