三维离散纵标-蒙特卡罗耦合系统TDOMINO开发与验证

韩静茹，陈义学，袁龙军，张春明

(1.华北电力大学 核科学与工程学院，北京 102206；2.环境保护部 核与辐射安全中心，北京 100082)

在求解辐射屏蔽问题时，常用的蒙特卡罗方法和离散纵标法各有长处和局限性。三维离散纵标(SN)-蒙特卡罗(MC)耦合是利用粒子角通量密度分布和MC源变量抽样概率之间的转换算法，将三维离散纵标程序与蒙特卡罗程序耦合，充分发挥各自的优势，同时又克服各自的局限性。对大型复杂屏蔽问题而言，三维耦合计算可更加有效地解决深穿透问题和详细描述复杂几何，揭示更为真实的安全裕量，有利于核电站安全性的改进和经济性的提高。自20世纪70年代至今，国内外科研人员对离散纵标-蒙特卡罗耦合计算进行了大量的研究，如美国橡树岭国家实验室开发的DOMINO[1]、美国加利福尼亚大学开发的PROBGEN[2]、国内西安交通大学开发的DORT2MCNP[3]等，但上述研究仅支持r-z几何下的二维离散纵标计算和蒙特卡罗计算的耦合。日本东芝公司实现了三维离散纵标程序TORT和MCNP的耦合[4]，但仅限于x-y-z几何，未解决r-θ-z几何问题。这严重限制了耦合方法在实际工程中的应用。因此，本文进行三维离散纵标-蒙特卡罗耦合系统TDOMINO的开发，并针对国际通用基准题进行验证计算。

1 三维SN-MC耦合系统框架

TDOMINO耦合系统流程图示于图1。以现有程序为基础，通过开发建立相应的耦合接口和源粒子抽样程序来实现程序数据交换和计算。这样的耦合框架既避免了对现有程序结构的修改，也便于今后对程序进行维护和升级。同时也使耦合形式具有较好的灵活性，可根据问题分析的需要选择用于耦合计算的程序。构成TDOMINO耦合系统的基础程序包括：三维离散纵标程序TORT[5]和蒙特卡罗程序MCNP[6]。在耦合系统中，为充分发挥两种程序的优势，TORT和MCNP分别用以模拟大块屏蔽区和复杂几何区计算。其中三维SN计算得到连接面的粒子角通量密度，并以BNDRYS格式存储输出，作为SN-MC接口程序的输入文件之一。SN-MC接口程序将连接面的角通量密度转换为粒子概率分布，再经MC-SOURCE源粒子抽样程序转换为MC源粒子信息参数，为MC计算提供源项，实现三维SN和MC耦合输运计算，得到复杂几何区的粒子通量分布。该程序系统可处理三维x-y-z及r-θ-z几何。

图1 TDOMINO耦合系统流程图

2 三维SN-MC耦合方法

对于上述两个基础程序，由于其物理建模和数值求解方法不同，耦合需解决连接区内SN计算的粒子角通量密度分布和MC计算所需源项之间映射的问题。在TDOMINO耦合系统中，通过开发三维SN-MC接口程序和MC-SOURCE源粒子抽样程序分别实现连接区内SN计算得到的粒子角通量密度和粒子概率分布之间映射及粒子概率分布和MC源粒子信息间的转换。具体算法为：

(1)

(2)

(3)

(4)

式中：ψ(ri，Eg，Ωm)为与空间、能量和方向相关的粒子角通量密度；i为空间网格(总网格数用I表示)；g为能群 (总能群数用G表示)；m为离散方向(总离散方向数用M表示)；Ai为第i个空间网格的面积；λm为SN求积组的方向余弦μm、ξm或ηm；wm为λm对应的SN求积组的权重；p(i)、p(g/i)和p(m/g/i)分别为粒子落在空间网格区间ri内的概率、网格区间ri内粒子能量落在能群Eg内的概率和网格区间ri、能群Eg内粒子落在离散方向m内的概率。

MC-SOURCE源粒子抽样程序基于三维SN-MC接口程序计算得到的粒子概率分布，根据式(4)进行随机抽样，依次获得粒子所在的SN计算划分的网格Δr、能群ΔE和离散方向ΔΩ区间，并在这些子区内均匀抽样，最终依次确定源粒子的位置(x,y,z)、能量(E)及飞行方向与(x,y,z)3个坐标轴的夹角余弦值(u，v，w)等信息，为下一步MC计算提供源项。需注意，圆柱坐标系下的源粒子位置和飞行方向表达方式与直角坐标系下的不同，需进行相应的坐标转换。

3 耦合程序的验证

为验证三维SN-MC耦合程序，采用美国橡树岭国家实验室公布的HBR-2屏蔽基准题[7]，利用TDOMINO耦合系统进行验证计算，并针对圆柱坐标系和直角坐标系下的耦合模型特点，采用不同的耦合形式进行计算。

三维SN-MC耦合计算可分为以下步骤。1) 模型划分：分为适合SN计算的屏蔽区及适合MC模拟的复杂几何区；2) 指定公共几何:连接大块屏蔽区(SN网格区)和复杂几何区(MC模型)的公共面；3) SN计算：得到公共面内每个SN网格的角粒子通量密度分布；4) 转换计算：将SN网格的角通量密度分别转换为粒子在空间、能量和方向上的分布；5) MC源抽样概率：将粒子分布转换成MC源参数相应的概率分布，并嵌入MC源抽样程序；6) MC计算：根据上步得到的源参数分布概率进行抽样模拟，进行MC计算。

3.1 HBR-2基准题

HBR-2基准题是在美国橡树岭国家实验室发布的用于验证压水堆压力容器屏蔽计算例题。HBR-2是热功率为2 300 MW的压水堆核电机组，反应堆堆芯包含157个燃料组件，压力容器内半径为197.485 cm。堆芯外依次围绕有堆芯围板、吊篮、热屏蔽、辐照监督管、压力容器以及生物屏蔽等部件，更详细的描述见HBR-2基准报告[7]。辐照监督管中心与堆芯中心的距离为191.15 cm，设置在热屏蔽外壁，位于20°方位角处。基准实验在反应堆第9循环内辐照监督管处放置测量仪对试样比活度进行测量。基准报告中提供了实验测量值和1/8模型的DORT程序计算结果。

3.2 计算模型

根据基准题提供的具体参数及反应堆对称性，结合TDOMINO应用特点，分别建立r-θ-z几何模型和x-y-z几何模型进行耦合计算。堆芯部分采用三维SN程序TORT建模计算，热屏蔽外的辐照监督管采用MCNP程序建立精确模型。

r-θ-z坐标下三维SN-MC耦合计算模型示于图2。以基准报告中计算模型为例，建立1/8堆型模型，0°和45°的角边界设为反射边界，压力容器外表面设为真空边界，同时将堆芯上反射层和下反射层的外表面设为轴向真空边界条件。将吊篮外表面设为SN-MC耦合模型的连接面，考虑到中子的散射作用，SN模型和MC模型有一段重叠区。堆芯到热屏蔽设为SN模型，采用TORT圆柱坐标系建模计算，吊篮外表面到压力容器设为MC模型，采用MCNP程序详细描述辐照监督管模型及材料，并在辐照监督管处设置了中子注量率计数卡。TDOMINO计算中，TORT和MCNP程序分别采用基于ENDF/B-Ⅵ的TEXT-10多群截面库和点截面库。耦合计算中TORT在r、θ、z三个方向上的网格个数分别取为118、52、88。

图2 三维SN-MC耦合计算r-θ-z模型(1/8堆芯)

为验证TDOMINO在x-y-z几何中的应用及处理多面耦合的能力，结合TORT建模特点，建立的1/4 HBR-2全堆芯组件均匀模型为三维SN-MC耦合计算模型，如图3所示。将围板设为耦合面，堆芯到围板的SN模型区采用TORT程序建立x-y-z几何模型(图3b)，围板到压力容器采用MCNP程序建模计算(图3c)，实现直角坐标系下的三维SN-MC多面耦合计算。

a——耦合计算模型；b——TORT计算模型；c——MCNP计算模型

3.3 计算结果

分别采用三维SN-MC耦合程序TDOMINO和单独的MCNP程序对HBR-2基准模型辐照监督管处6种典型核素的比活度进行计算，并与实验测量、离散纵标法程序计算结果进行对比分析。

表1列出HBR-2基准模型辐照监督管处6种典型核素比活度的实验测量值。表2列出不同程序计算的6种核素比活度与实验测量值间的比值(C/M)。其中实验测量值和DORT采用BUGLE-96[8]库的计算结果取自HBR-2基准报告，TORT结合TEXT10多群截面数据库的计算结果取自文献[9]，SN-MC和SN-MC-2分别表示TDOMINO在圆柱坐标系和直角坐标系下的计算结果。MCNP表示堆芯组件均匀模型的单一MC计算结果，MCNP-PIN表示堆芯PIN功率的计算结果取自文献[10]。DORT、SN-MC、SN-MC-2、MCNP、MCNP-PIN和TORT程序计算得到的辐照监督管处的平均C/M值分别为0.89±0.04、1.03±0.04、1.04±0.03、1.10±0.04、0.91±0.03和1.04±0.04。计算结果表明，TDOMINO耦合系统计算结果与实验测量值和其他程序计算值吻合较好，初步验证了三维SN-MC耦合系统TDOMINO在不同空间坐标系中应用的正确性。在上述计算中，SN-MC耦合程序的计算时间为3 908 min，其中，SN-MC中MCNP计算耗时3 869 min，最大统计误差为3%左右，SN-MC中TORT计算耗时39 min。单一MCNP程序计算耗时7 180 min，最大统计误差为5%左右。表明三维SN-MC耦合程序与单一MCNP程序相比，在较少的时间下得到的统计误差更小，计算结果可靠性更强。

表1 实验测量所得的比活度

表2 计算与测量的比活度之比

4 结论

为校核验证三维离散纵标-蒙特卡罗耦合程序系统在复杂模型辐射屏蔽计算问题中的应用可行性，针对美国HBR-2压力容器基准实验，对辐照监督管处的中子能谱和几种典型核素的比活度进行了计算分析。采用三维TDOMINO耦合系统分别建立了圆柱坐标系下的HBR-2 1/8模型和直角坐标系下的HBR-2 1/4耦合计算模型，并进行了验证计算。在辐照监督管处计算的平均C/M值分别为1.03±0.04和1.04±0.03，可看出耦合计算结果与实验测量及其他程序计算值吻合良好，证明了TDOMINO耦合系统开发的正确性。下一步将针对具体的工程实际问题(如反应堆堆腔漏束辐射计算等)，开展三维离散纵标-蒙特卡罗耦合系统的应用研究。

参考文献：

[1] EMMETT M B, BURGART C E, HOFFMAN T J. DOMINO, a general purpose code for coupling discrete ordinates and Monte Carlo radiation transport calculations, ORNL-4853[R]. USA: ORNL, 1973.

[2] EGGLESTON J E, ABDOU M A, YOUSSEF M Z. The use of MCNP for neutronics calculations within large buildings of fusion facilities[J]. Fusion Engineering and Design, 1998, 42(1-4): 281-288.

[3] 郑征，吴宏春，曹良志，等. 蒙特卡罗与离散纵标耦合方法在压水堆堆腔漏束计算中的应用[J]. 强激光与粒子束,2012,24(12):2 946-2 950.

ZHENG Zheng, WU Hongchun, CAO Liangzhi,et al. Application of Monte Carlo and discrete ordinate coupling method to pressurized water reactor cavity radiation streaming calculation[J]. High Power Laser and Particle Beams, 2012, 24(12): 2 946-2 950(in Chinese).

[4] KUROSAWA M. TORT/MCNP coupling method for the calculation of neutron flux around a core of BWR[J]. Radiation Protection Dosimetry, 2005, 116(1-4): 513-517.

[5] BRIESMEISTER J F. MCNP: A general Monte CarloN-particle transport code, Version 4C, LA-13709-M[R]. USA: LANL, 2000.

[6] RHOADES W A, SIMPSON D B. The TORT three-dimensional discrete ordinates neutron/photon transport code, ORNL/TM13221[R]. USA: ORNL, 1997.

[7] REMEC F B K, KAM H B. Robinson-2 pressure vessel benchmark[R]. USA: ORNL, 1997.

[8] WHITE J E. BUGLE-96: Coupled 47 neutron, 20 gamma-ray group cross section library derived from ENMF/B-Ⅵ for LWR shielding and pressure vessel dosimetry applications[R]. USA: RSIC Data Library Collection, 1996.

[9] 杨寿海,陈义学,王伟金,等. 三维离散纵标方法在RPV快中子注量率计算中的初步应用[J]. 核科学与工程，2011，31(4)：294-298.

YANG Shouhai, CHEN Yixue, WANG Weijin, et al. The analysis of RPV fast neutron flux calculation for PWR with three-dimensional SNmethod[J]. Chin J Nucl Sci Eng, 2011, 31(4): 294-298(in Chinese).

[10] 石生春. 基于蒙特卡罗方法的压水堆压力容器快中子注量率的计算分析[D]. 北京：华北电力大学，2010.