研究堆退役中的快速辐射传输计算方法研究

毕远杰，杨宏伟，郭金森，骆志平，李传龙，陈 凌

（中国原子能科学研究院 辐射安全研究所，北京 102413）

研究堆退役是核燃料循环的重要环节之一，目前我国的研究堆已逐步进入退役阶段。然而，现有的退役工作通常是基于初步的源项分析和辐射监测数据进行的。这样的操作流程存在的问题是：随着设备的拆除和去污工作的进行，辐射源项和屏蔽条件均在不停地发生变化，而辐射监测数据不能及时地跟进。同时，在工作过程中，即使选取了参考点进行辐射剂量的测量，由于辐射场分布的差别，参考点也不一定能反映人员所在位置的实际剂量。因此，操作人员可能在实际的工作过程中受到额外的照射，对他们的身体造成伤害。

随着计算机技术的发展，目前在大型核设施的退役过程中，人们更希望能通过计算机模拟给出整个厂房内部的三维剂量场分布及其动态变化情况，为工作人员在放射性工作场所内进行放射性操作时的工作方案设计提供保障。根据计算机的模拟结果，工作人员可对不同的退役方案进行评估，选择最优化的方案进行实施。因而，进行三维辐射场快速计算成为时代的需求。

在进行三维辐射场快速计算的过程中存在两个问题：1）空间，需进行整个厂房内全空间的辐射场计算；2）时间，伴随着设备的拆卸需在尽量短的时间内给出新的剂量分布情况。目前传统的计算方法很难同时解决这两个问题，因此，本文拟制定一种新的计算方案，可在较短的时间内给出全空间的辐射场分布情况。

1 辐射场计算方法现状

1.1 传统计算方案

解决辐射传输问题的基本工具是Boltzmann方程［1-2］。求解Boltzmann 方程的方法分为两类：1）确定论方法，目前最常用的是离散纵标方法，代 表 性 的 程 序 有ANISN［3］、DORT［4］、TORT［5］等；2）蒙特卡罗方法（简称蒙卡方法），蒙卡方法是一种随机性方法，代表性的程序有MCNP［6］、FLUKA［7-8］、GEANT4［9-10］、MARS［11］、PHITS［12］等。在实际的工程应用中，还常用到Boltzmann方程的近似解法，如在屏蔽计算中用点核方法进行估算，常用的点核程序有MICROSHIELD［13］、QAD［14］、MERCURE［15］等。

点核方法的主要优点是快速。除广泛用于屏蔽设计的估算外，近年来新兴的用于辐射防护最优化（ALARA）的辐射场快速计算软件也多基于点核方法编写。国际上已有一些相关软件，如比利时核研究中心开发的VISIPLAN 软件［16］；法 国 的CEA-LIST 在 点 核 程 序MERCURE 5的基础上开发的NARVEOS软件［17］；日本核燃料循环开发署与日本原子能研究机构、挪威OECD／NEA Halden的反应堆项目联合开发的退役工程支持系统DEXUS［18］；美国电力研究协会EPRI 开发的三维可视化系统［19］等。国内方面如中国科学院核能安全技术研究所FDS团队基于国家大科学工程全超导托卡马克核聚变实验装置EAST 建立的辐射虚拟仿真系统RVIS［20］，也将建立基于点核方法的快速计算模块。以上基于点核方法的虚拟仿真软件虽快速，但存在源项周围剂量估算不准确、不能处理复杂的几何结构、可靠性低等问题。

1.2 辐射场耦合计算

蒙卡、离散纵标和点核方法各有优缺点（表1）。单一的计算方法均存在各自的不足，蒙卡方法对于大空间厚屏蔽问题无法在合理的时间内给出可靠的计算结果，确定论方法则不适合复杂源项和几何区域。耦合计算则可将整个问题划分为两个区域进行求解，在复杂源项和几何区域采用蒙卡方法进行模拟，在大空间厚屏蔽区域采用确定论方法求解深穿透问题。

1）蒙卡-离散纵标耦合方法

目前存在的蒙卡-离散纵标耦合方法的研究，主要是针对MCNP和DOORS代码包中程序的耦合。如MCNP 与一维离散纵标程序ANISN 的 耦 合［21-23］，MCNP 与 二 维 离 散 纵 标程 序DORT 的 耦 合［24-25］，MCNP 与 三 维 离 散纵标程序TORT 的耦合［26-27］。蒙卡-离散纵标耦合的特色在于进行全空间的趋于精确的计算。但在离散纵标方法中，考虑到能量多群、方位角离散和空间网格划分，得到的是一系列高维积分微分方程组，需各种技巧迭代至收敛。尤其是三维离散纵标程序，计算量庞大。同时，在蒙卡和离散纵标的耦合过程中，若连接面上通量的空间分布和方位角分布缺乏对称性，则需调整方位角离散和空间网格划分的大小来保证离散纵标方法求解的精度，这会大幅增加计算量。所以，蒙卡-离散纵标耦合方法并不适用于快速计算。

表1 3种计算方法的优缺点Table 1 Comparison of advantage and disadvantage of three methods

2）蒙卡-点核耦合方法

蒙卡与点核耦合进行全空间辐射场快速计算还未见详细报道，蒙卡-点核耦合方法的特点为：（1）可给出全空间的三维辐射场，解决点核方法在复杂源项和几何区域存在的问题；（2）相对于蒙卡-离散纵标耦合方法的趋于精确计算，蒙卡-点核耦合方法更加快速，可随着设备的拆卸在尽量短的时间内给出新的剂量场分布。

本文基于蒙卡-点核耦合的计算方案，编写耦合接口模块，用于辐射场的三维快速计算。

2 耦合计算系统

图1为耦合计算程序系统示意图。在耦合计算中，将整个问题划分为两个区域：复杂源项和几何区域及简单几何大屏蔽区域。复杂源项和几何区域的结构通常较复杂，所以在这种区域采用蒙卡方法进行粒子输运模拟。而在大屏蔽部分，几何通常较简单，同时属于粒子深穿透问题，采用点核方法进行计算。结合蒙卡程序和点核计算的结果，得到三维辐射场快速计算软件，并对软件进行基准校验和对比分析。

图1 耦合计算程序系统示意图Fig.1 Schematic of coupled calculation program

耦合接口模块的编制是计算程序系统的重点，耦合接口模块需解决的主要问题是如何将蒙卡输出的粒子轨迹转换为点核计算的源项。首先需从物理上建立两者之间的联系，即获得等效面源。通过分析Boltzmann方程可得出：在一个封闭曲面内的源和物质，可由这个封闭曲面上的粒子流量率来代替，这两种情况对探测器的响应相同［1］。在这个等效问题中，通过连接曲面的流量率相当于等效面源的源强。

在计算过程中，首先对交界面进行空间上的网格划分得到面元dA，接着对dA 出射的粒子进行方位角离散和能量分群，最终得到（dA，dE，dΩ）。设dA 的中心点为A′、外法向为nA，则位于P 点的探测器和dA 的外法线方向的夹角，因 此，dA 对P点探测器的贡献为：

其中：d（θ）为沿θ方向穿过屏蔽的厚度；r为点到探测器的距离；Hθ为无屏蔽的情况下在θ 方向与源单位距离处的剂量率；μ 为γ射线的线衰减系数；B（d）为累积因子。

3 退役算例

本文对退役过程中放射性废物储存容器所在的大尺度厚屏蔽房间的辐射场进行了计算，计算算例几何和参考点如图2所示。圆柱形放射性废物储存罐位于x＝-290cm 处，内部装有含放射性核素137Cs的淤泥，容器壁为1cm厚不锈钢。放射性核素浓度在z方向呈指数分布，抽样函数为：其中：淤泥在z 方向的范围为从z1到z2；ξ 为［0，1］上均匀分布的随机变量；q为系数。

图2 计算算例几何和参考点Fig.2 Geometry of model and reference points

在x＝-683cm～x＝-663cm 放置一短混凝土墙，在x＝-1 283cm～x＝-1 263cm放置一长混凝土墙。

计算结果采用通用蒙卡软件FLUKA 进行校准，表2比较了FLUKA 和耦合程序的计算结果。对于选择的参考点，两个计算程序的相对偏差小于10%。在计算时间方面，将整个房间划分为105个网格进行三维辐射场计算，在3.4GHz、4GB内存的个人计算机上进行测试，FLUKA 达到表2 中的计算误差所需的计算时间为12h，耦合程序用20min得到相对误差1%的等效面源，18min进行点核计算。

表2 参考点周围剂量当量率比对Table 2 Comparison of ambient dose equivalent in reference points

4 结论

本文基于蒙卡与点核耦合的辐射场计算方法，编制了相应的耦合计算程序，用于研究堆退役过程中的辐射场快速计算。蒙卡-点核耦合方法结合了蒙卡和点核计算各自的优点，不仅适合于复杂源项和屏蔽体区域的计算，同时也适合于大屏蔽区域的辐射场快速计算。作为算例，本文采用编制的耦合程序对研究堆退役过程中放射性废物储存房间的辐射场进行了计算，计算结果验证了程序的可靠性和适用性。

［1］ CHILTON B，SHULTIS J K，FAW R E.Principles of radiation shielding［M］.USA：Prentice-Hall，Inc.，1984.

［2］ SHULTIS J K，FAW R E.Radiation shielding［M］.USA：American Nuclear Society，Inc.，2000.

［3］ ENGLE W W，Jr.ANISN，a one-dimensional discrete ordinates transport code with anisotropic scattering，K-1693［R］.USA：ORNL，1967.

［4］ RHOADES W A，CHILDS R L.The DORT two-dimensional discrete ordinates transport code［J］.Nuclear Science and Engineering，1988，99：88-89.

［5］ AMY Y Y.The three-dimensional，discrete ordinates neutral particle transport code TORT：An overview［C］∥OECD／NEA Meeting on 3D Deterministic Radiation Transport Computer Programs，Feature，Applications and Perspectives.Paris：［s.n.］，1996.

［6］ HENDRICKS J S，MCKINNEY G W，WATERS L S，et al.MCNPX user's manual：Ver-sion 2.5.0，LAUR-05-2675［R］.USA：LANL，2005.

［7］ BOHLEN T T，CERUTTI F，CHIN M P W，et al.The FLUKA code：Developments and challenges for high energy and medical applications［J］.Nuclear Data Sheets，2014，120：211-214.

［8］ FASSO A，FERRARI A，RANFT J，et al.FLUKA： A multi-particle transport code，SLAC-R-773［R］.［S.l.］：［s.n.］，2005.

［9］ AGOSTINELLI S，ALLISON J，AMAKO K，et al.GEANT4：A simulation toolkit［J］.Nucl Instrum Phys Methods Res A，2003，506（3）：250-303.

［10］ALLISON J，AMAKO K，APOSTOLAKIS J，et al.GEANT4 developments and applications［J］.IEEE Trans Nucl Sci，2006，53（1）：270.

［11］MOKHOV N V，STRIGANOV S I.MARS15 overview［C］.［S.l.］：［s.n.］，2007.

［12］NIITA K，MATSUDA N，IWAMOTO Y，et al.PHITS overview［C］.［S.l.］：［s.n.］，2007.

［13］MICROSHIELD®user's manual［M］.［S.l.］：Grove Software，INC.，2009.

［14］CAIN V R.QAD-CG，a combinatorial geometry version of QAD-P5A，apoint kernel code for neutron and gamma-ray shielding calculations，ORNL-CCC307［R］.USA：ORNL，1977.

［15］DUPONT C，NIMAL J C.MERCURE-4：A three-dimensional Monte Carlo program for the integration of specific point attenuation kernels in a straight line，OLS-82-116［R］.［S.l.］：［s.n.］，1980.

［16］VERMEERSCH F.Dose assessment and dose optimization in decommissioning using the VISIPLAN 3D ALARA planning tool［C］∥Radiation Protection and Decommissioning Conference ABR／BVS.Brussels：［s.n.］，2003.

［17］THEVENON J B.Narveos：A new tool supporting ALARA studies［R］.French：Marcoule Cea Center，2011.

［18］IGUCHI Y，KANEHIRA Y，TACHIBANA M，et al.Development of decommissioning engineer-ing support system （DEXUS）of the Fugen Nuclear Power Station［J］.J Nucl Sci Technol，2004，41（3）：367-375.

［19］Work planning and 3D visualization tool［R］.USA：Electric Power Research Institute，2009.

［20］TANG Z D，LONG P C，HUANG S Q，et al.Real-time dose assessment and visualization of radiation field for EAST Tokamak［J］.Fusion Engineering and Design，2010，85（7-9）：1 591-1 599.

［21］GALLMEIER F X，PEVEY R E.Creation of a set of interface utilities to allow coupled Monte Carlo／discrete ordinate shielding analysis［C］∥The Third International Topical Meeting on Nuclear Applications of Accelerator Technology.USA：American Nuclear Society，1999：404-409.

［22］GABRIEL T A.CALOR：A Monte Carlo program package for the design and analysis of calorimeter systems，ORNL／TM-5619［R］.USA：ORNL，1977.

［23］CLOTH P，FILGES D，NEEF R D，et al.HERMES，a Monte Carlo program system for beam-materials interaction studies［R］.Germany：KFA，1988.

［24］JOHNSON J O.Shielding design of the spallation neutron source，shielding aspects of accelerators，targets and irradiation facilities［C］∥SATIF 4 Workshop Proceedings.Knoxville：［s.n.］，1998：89-100.

［25］JOHNSON J O，SANTORO R T，LILLIE R A，et al.The SNS target station preliminary title I shielding analyses［J］.J Nucl Sci Technol，2000，37（Suppl.）：35-39.

［26］CHEN Y.Coupled Monte Carlo discrete ordinates computational scheme for three-dimensional shielding calculations of large and complex nuclear facilities［D］.Forschungszentrum Karlsruhe GmbH：Institute Für Reaktorsicherheit，2005.

［27］韩静茹.三维蒙特卡罗-离散纵标双向耦合屏蔽计算方法研究［D］.北京：华北电力大学，2012.