韩静茹1.，陈义学1，袁龙军1
（1.华北电力大学核科学与工程学院，北京102206：2.环境保护部核与辐射安全中心，北京100082）
摘要：蒙特卡罗方法（MC）和离敏额标方续（SN）是辐射屏蔽设计靠用的两种计算方法.蒙特卡罗方法可以精确描述复杂模型，但对于深穿透问题计算耗时：离敏纵标方法计算速度较快，但在儿何描述方面存在限制，多维计算还存在射线效应，单一的MC或SN方法在大型复杂核装置屏藏计算中都存在一通过自主开发的接口程序实现MC粒子概率分布与SN角通量密度之间的相互转换，实现蒙特卡罗和 定限制，为了有效解决大需复杂核设链屏蔽计算问题，本文研究了蒙特卡罗-离散银标双向橘合方法.高散纵标方法的耦合计算，充分发挥两种方法的优势，本文进行了直角坐标系和医柱坐标系下的测试验证，耦合程序计算结果与MCNP、TORT结集吻合良好，初步证明了蒙特卡罗-离散织标双向耦合程序 的正确性.
关健词：蒙特卡罗：离散纵标：双向合
中图分类号：TL329 文繁标志码：A
文章编号：0258-0918（2013)04-0368-06
Developmentof Three-dimensionalProgramBasedon Monte Carlo and Discrete OrdinatesBidirectional CouplingMethod
HAN Jing-ru' CHEN Yi-xue' YUAN Long-jun
(1. School of Nucdlear Scienee and Engineering Noeth China Electric Power University Beijing 102206;2. Nuelear asd Radiation Sefety Center of MEP Bejing 100082)
Abstract; The Monte Carlo (MC) and discrete ordinates (SN) are the monly used methods in the design of radiation shielding. Monte Carlo method is able to treat thegeometry exactly but time-consuming in dealing with the deep penetrationproblem. The discrete ordinate method has great cotmputational efficiency but it is quitecostly in puter memory and it suffers from ray effect. Single discrete ordinates method or single Monte Carlo method has limitation in shielding calculation for large
ordinates bidireetional coupling method is developed. The bidirectional coupling method plex nuclear facilities. In order to solve the problem the Monte Carlo and discreteis implemented in the interface program to transfer the particle probability distributionof MC and angular flux of discrete ordinates. The coupling method bines thebeen calculated by the coupling methods. The calculation results are performed with advantages of MC and SN. The test problems of cartesian and cylindrical coordinate haveparison to MCNP and TORT and satisfactory agreements are obtained. Thecorrectness of the program is proved.
Key words: Monte Carlo; discrete ordinates; bidirectional coupling
辐射屏藏系统设计是核工程设计的核心内上存在限制，而且同样为单向耦合，为了有效究了蒙特卡罗-离散纵标双向耦合方法，并进行 了直角坐标系和画柱坐标系下的测试验证，同MCNP和TORT计算结果进行比较分析，初步验证了三维蒙特卡罗-离散纵标双向耦合方
容之一，其设计优劣直接影响核系统寿命及工求解大型三维复杂核装置的屏蔽计算问题，本作人员与周围环境的辐射安全.选取合适的屏文以蒙特卡罗方法和离散纵标方法为基础，研蔽计算方法是保证疆射屏蔽系统设计质量的关 键.蒙特卡罗方法（MC）和离散纵标方法（SN）是核工程设计中常用的屏蔽计算方法.MC方法的优点是几何描述精确，但计算耗时，尤其难以精确解决屏藏计算中常见的深穿透问题.现法的正确性和可行性.统计误差的连续传递性，使其在研究和实际应 有的MC分段计算方法由于分段计算中MC用中存在很大的困难.与MC方法相比，SN应.对于同时具有深穿透和复杂几何的屏蔽计子输运问题的两种不同方法，MC方法通过 方法计算速度快，但耗费内存，而且存在射线效算，单一的MC或SN方法难以提供可靠的计算结果.
1蒙特卡罗-离散纵标双向耦合方 法及程序
蒙特卡罗方法和离散纵标方法是求解粒模拟，然后用统计方法得出随机变量某个特
对大量单个粒子的运动历史逐个进行跟踪征的估计量，作为问题的解.SN方法是在相空间（r×EX0）内对空间r、能量E和方向变量Ω采用直接离散方法数值求解中子 输运方程，得到网格△r、能群△E和离散方向△区间内的粒子角通量密度y（r，E，Ω）.为了克服两种计算方法在实际工程应用中的困难，结合蒙特卡罗方法和离散纵标方法的优势，开发三维鹏合方法，实现MC 和SN方法之间的耦合计算，关键在于蒙特卡罗方法模拟的粒子径迹信息和离散纵标方法计算的角通量密度之间的相互转换.其对应转换方法基于自主研发的三维映射 方法[1，如式（1)-（4）所示.
基于MC和SN方法的耦合方法可以有效解决同时具有深穿透和复杂几何的屏蔽问题，基于 这种思想，国内外开展了大量的研究，取得了一定的研究成果，如耦合程序MCNP-ANISN）、MCNP-DORT2、HETC96-ANISN、HERMESANISN、MCNP-TRIDENT、DORT- PROBEGEN-MCNPC.DORT-DOMINO-MORSEDORT-MCNP等.然而这些合方法均为单向耦合，仅支持MC方法与一维或二维SN耦合.这些特点严重限制了合方法在实际屏蔽计算中的应用.此外， 陈义学等开发了三维McDeLicious-TORT耦合程序，并用于直角坐标系下国际聚变材料辐照装置（IFMIF）屏蔽设计计算.日本东芝公司的MasahikoKurosawa实现了直角坐标系 下的TORT-MCNP1耦合.上述三维MCSN与SN-MC耦合程序在圆柱坐标系的应用
d8（E E.）8（，F.）8（0.a.）（1)
方向.式（2）~式（4）表示SN角通量密度到MC粒子径信息的转换关系.其中，（r，E）表示网格能群E和离散方向区 间内的粒子角通量密度；入表示粒子飞行方向与面法线方向夹角的余弦值：A表示网格区间对应的连接面的面积；求和上限I、G和M分别表示连接面网格数目、能群数目和离散方向数目：p（r）表示MC抽样粒子位于网格区间内的概率；p（E/）表示在网格区间内抽样，粒子位于能群区间E内的概率；p（/E/r）表示在网格区间，能群区间E内抽样，粒子位于离散方向0内的概率.
(2)
(3)
(4)
式（1）表示MC粒子径迹信息到SN角通量密度的转换关系.其中，表示空间网格（i.j.k）中，第g能群内，m方向范围内的粒子角通量密度：ueight表示MC粒子n的权 重；w表示求积权重系数：N表示MC模拟的源粒子数；△A表示给定面元的面积；A.表示MC粒子径迹和面法线方向夹角的余弦值：E.，和0.分别表示粒子的能量、空间位置和飞行
三维蒙特卡罗-离散纵标双向耦合计算程序系统流程图如图1所示.其中蒙特卡罗模拟采用国际通用的粒子输运程序MCNP，离散纵标计算采用美国橡树岭国家实验室开发的三维离散纵标法程序TORT，面接口程序则基 于上述三维映射方法，实现MCNP模拟的粒子径迹信息与TORT计算得到的粒子通量分布之间的相互转化.
图1三维蒙特卡罗-离散纵标双向耦合计算程序系统流程图
Fig. 1 Flow chart of the program system for three-dimensional bidirectional coupled MC-SN
为44×20×20.MCNP、TORT、MC-SN、SN-MCMCNP结果相比，其他四种程序计算结果相差不 和MC-SN-MC五种程序计算结果见表1.与超过2%.其中MCNP计算时间为793分钟，计算结果的统计误差小于0.5%.相近统计误差下，MC-SN-MC计算时间为427分钟.X-10cm处 探测器的中子能谱如图3所示.通过对比可以看出，计算结果吻合良好，初步验证了直角坐标系下蒙特卡罗-离散纵标双向耦合程序系统的正确性和计算中选取不同耦合方案的可行性.
2三维耦合程序系统测试计算
为了验证三维蒙特卡罗-离散纵标双向耦合程序系统的正确性和可行性，选取直角坐标系和圆柱坐标系下的两个简单模型，分别采用MC-SN.SN-MC及MC-SN-MC三种不同的耦 合方案对其进行屏蔽计算，并与单独的MCNP计算和TORT计算结果进行对比.三维耦合程序系统中，蒙特卡罗方法计算采用MCNP程序及点截面数据库，三维离散纵标方法计算采 用TORT程序及TEXT-10多群截面数据库.
2.1X-Y-Z模型测试
直角坐标系下的测试模型为长20cm、宽20cm、高25cm的长方体，其中底端为20cm×20cm×6cm的中子源，探测器位于模型上端的 计数区，沿X轴方向共设置20个探测器进行计数.为了充分发挥SN方法快速解决深穿透问题的优势，MC-SN和SN-MC耦合面分别设为距底端10cm和18cm处的水平面，如图2所 示.中子源为均匀分布各向同性的14MeV中子，屏蔽材料为不锈钢和水的混合物（60a%SS31640a%Water).
图2X-Y-Z测试模型
TORT计算P3阶勒让德展开，S8全对称高斯求积组，离散方向个数为96，网格划分数目
Fig.2Test model of X-Y-Z
表1不同程序计算得到的探测器处中子注量率分布
Table 1Neutron flux distrilbution in the detects according to different programs
x/cm MCNP TORT MC-SN SN-MC MC-SN-MC1 1. 195 7E-03 1. 192 8E-03 1. 190 7E-03 1. 195 5E-03 1.194 3E-032 1. 193 8E-03 1.192 8E-03 1. 184 0E-03 1. 192 0E-03 1.193 8E-033 1.198 2E-03 1. 190 5SE-03 1. 192 8E-03 1. 192 8E-03 1. 186 6E-03 1. 192 0E-03 1. 192 1E-03 1. 175 7E-03 1. 189 7E-03 1.193 4E-035 4 1. 197 1E-03 1. 192 8E-03 1. 194 4E-03 1. 188 3E-03 1.198 4E-036 1. 199 8E-03 20-382611 1. 192 0E-03 1.196 2E-03 1.194 4E-037 1.191 SE-03 1 199 0E-03 1. 192 8E-03 1. 184 4E-03 1. 193 6E-03 1.190 9E-038 9 1. 195 1E-03 1.192 8E-03 1. 192 8E-03 1. 190 0E-03 1. 189 5E-03 1. 186 5E-03 1.195 6E-03 1. 191 3E-03 1.190E-0310 1. 196 2E-03 1.192 8E-03 1. 192 7E-03 1. 200 9E-03 1.196 SE-0311 1.200 5E-03 1.192 8E-03 1.194 3E-03 1. 185 0E-03 1.196 5E-0312 13 1. 194 6E-03 1. 193 9E-03 1. 192 8E-03 1. 192 8E-03 1.189 0E-03 1.187 0E-03 1. 198 3E-03 1.188 9E-03 1.196 4E-03 1.189 7E-0314 1.187 9E-03 1.192 8E-03 1. 195 7E-03 1. 192 6E-03 1.186 7E-03
x/cm 15 1. 193 5E-03 MCNP 1.192 8E-03 TORT 1. 192 0E-03 MC-SN 1. 187 6E-03 SN-MC 1.195 0E-03 MC-SN-MC16 1. 196 9E-03 1. 192 8E-03 1. 195 2E-03 1. 198 0E-03 1. 189 7E-0317 1.193 3E-03 1.192 8E-03 1.193 8E-03 1.194 6E-03 1.189 5E-0318 1.190 7E-03 1.192 8E-03 1.190 3E-03 1.197 1E-03 1. 196 6E-0319 20 1.194 4E-03 1. 194 7E-03 1.192 8E-03 1. 192 9E-03 1.188 7E-03 1.187 7E-03 1.2048E-03 1. 193 1E-03 1. 192 2E-03 1. 191 8E-03
图4R-6Z测试模型Fig. 4Test model of R-0-Z
图3X-Y-Z模型中子能谱对比Fig.3Comparison of the neutronspectrum in X-Y-Z model
0.4%.相近统计误差下，MC-SN-MC计算时间为134分钟.角度为21.375处探测器的中算结果吻合良好.初步验证了圆柱坐标系下蒙 子能谱图如图5所示.通过对比可以看出，计特卡罗-离散纵标双向耦合程序系统的正确性和计算中选取不同耦合方案的可行性，
2.2R-0-Z模型测试
圆柱坐标系下的测试模型为半径为25cm、角度为45、高度为40cm的圆柱体，其 中有半径为6cm、角度为45°、高度为40cm的圆柱中子源，而计数区位于半径为21cm处轴向平面位置，沿周向设置20个探测器进行计数.MC-SN耦合面设为R=10cm处的平面，SN-MC耦合面设为R=18cm处的平面，如图 4所示.其中子源为均匀分布各向同性的14MeV中子.屏蔽材料为不锈钢和水的混合物(60a%SS31640a%Water).
对称高斯求积组，离散方向个数为96，网格划 TORT计算采用P3阶勒让德展开，S8全分数目为44×20×20.MCNP、TORT、MCSN.SN-MC和MC-SN-MC五种程序计算结果见表2.其他程序计算结果与MCNP计算结算时间为349分钟，计算结果的统计误差小于 果相比，最大误差不超过2%，其中，MCNP计
图5R-9-Z模型中子能谱对比Fig.5Comparison of the neutronspectrum in R-0-Z model
