郑 剑，李 斌，邹 俊，杨 琪

(中国科学院核能安全技术研究所，中国科学院中子输运理论与辐射安全重点实验室，合肥230031)

基于网格计数的停机剂量率计算方法研究

郑 剑，李 斌，邹 俊，杨 琪

(中国科学院核能安全技术研究所，中国科学院中子输运理论与辐射安全重点实验室，合肥230031)

本文根据严格二步法的计算理论，研究了基于网格计数的停机剂量率计算方法，设计并实现了基于网格计数的停机剂量率计算程序。该程序能够支持圆柱坐标下的网格计算。本文使用源子程序进行复杂源描述。为了加快计算速度，本文采用了多节点和多线程等技术。本文利用国际热核聚变实验堆(ITER)停堆剂量基准实验ITER-T426进行测试，计算结果与实验值吻合良好，证明了该方法和程序的正确性和可用性。

网格计数；停机剂量率；严格二步法；自动耦合

核装置运行停止或运行间歇期间，材料活化所释放的衰变光子对需要靠近或进入装置内部进行实验测量或维修检测工作的人员造成严重的辐射伤害。因此准确地评估核装置的停堆剂量率水平，对于核装置屏蔽设计、维修计划的制订以及核装置的退役有着重要的参考意义。

本文在FDS团队自主研发的大型集成多功能中子学计算分析系统VisualBUS[3]及已经开发的停机剂量率计算程序的基础上[4-5]，研究了基于网格计数的严格二步法，设计并实现了基于网格计数的停堆剂量率计算程序。该程序支持圆柱坐标下的网格计算，使用源子程序及并行技术，解决了传统停机剂量率计算过程中的精度不够、模型划分困难、源栅元数受限等问题，实现了核装置停机剂量率的精确计算，计算结果可以为装置的屏蔽设计、实验方案的制订以及工作人员允许靠近装置的等待时间等提供重要的参考依据。

1 计算方法

1.1 总体设计

网格计数停机剂量率计算方法基于栅元计数严格二步法的思想，充分利用MCNP程序的网格计数能力。图1描述了该方法的计算流程。

图1 计算流程图Fig.1 Calculation flow chart

该方法首先使用MCNP程序[6]，对已划分了网格的计算模型进行网格计数中子输运计算，获得所有被计数网格的中子分群通量。然后调用MCNP程序对每一个网格进行PTRAC计算，获得每个网格的粒子分布情况，用于计算网格的平均材料。

然后，该方法根据上面得到的网格内通量和材料，以及用户提供的辐照方案，通过使用FISPACT程序[7]对每一个计数网格进行活化计算，得到不同冷却时间下所有网格的衰变光子分布情况。之后，提取活化之后的光子谱及强度，生成光子源输入文件。

最后，使用MCNP程序对不同时间段下的衰变光子进行网格计数输运计算，最终获得不同时间段下计数网格的剂量率分布。

该方法通过接口程序，自动完成计算过程中信息的处理及传递，实现MCNP和FISPACT的双向自动耦合，避免了其中耗时的、复杂困难易出错的手工处理过程。

利用比浊法的原理，用光密度值(OD值)的变化来判断菌株存活情况。根据朗伯比尔定律OD=log It IO=kbc，c表示样品浓度，若样品液厚度b一定，则OD值与样品的浊度相关[11]。将活化好的菌种菌液以1%的接种比例，分别加入用1 mol/L HCl和1 mol/L NaOH调pH值(3，4，5，6，7)的灭菌YPD液体培养基中，28 ℃、180 r/min振荡培养48 h后，将菌液稀释至合适倍数(使吸光值在0.1～0.8之间)，在波长600 nm处测量吸光度值，每组试验做3组平行。

1.2 网格材料计算

MCNP程序提供了一种名为粒子径迹输出卡(PTRAC: Particle Track Output)的计算功能。通过在MCNP输入文件中填写PTRAC卡，并调用MCNP程序进行输运计算，得到名为ptrac的输出文件，该文件中打印出了计数粒子信息。统计网格中这些计数粒子所在的栅元信息，可以得知网格所含栅元情况，从而计算出网格的平均材料。计算公式如下：

(1)

由于PTRAC卡及源描述不支持如图2所示的圆柱几何划分出的扇形网格，本文研究并提出了一种方法来计算这类扇形几何网格内的材料：首先，该程序根据网格的边界设定了立方体型的包围体，在此包围体内进行PTRAC计算，然后去掉不在网格内的源粒子，接着依照上文所述方法计算网格材料组成。

图2 单元与包围体Fig.2 Cell and bounding volume

1.3 源子程序

本文开发了一种基于网格的衰变光子取样程序，并将该程序耦合到MCNP程序中，图3给出了源子程序抽样流程。

图3 源子程序流程图Fig.3 Source subroutine flowchart

源输入文件中包含了网格的几何坐标、体积、衰变光子强度、衰变光子谱、栅元号以及光子能群结构信息。源子程序使用随机数函数rang()在包围盒内进行抽样。

首先，源子程序按一定概率进行网格抽样。对于不同网格，其体积和光子强度值越大被抽样到的概率越大，抽样概率计算公式如下：

(2)

其中，Vi表示第i个网格的体积，Di表示第i个网格的光子强度，Ri表示从n个网格中抽取的第i个网格的体积与光子强度值。

然后，在抽样得到的网格内进行均匀光子抽样，通过公式(3)确定立方体网格内衰变光子的几何坐标(xxx,yyy,zzz)；通过公式(4)确定圆柱网格内衰变光子的半径r，然后通过坐标变换得到几何坐标(xxx,yyy,zzz)。对于得到的衰变光子坐标，调用MCNP程序提供的库函数确定该光子所在的栅元号icl，本文支持对重复结构的栅元处理。

x=(xmax-xmin)·rang()+xmin

(3)

(4)

最后，利用网格的光子谱和光子强度信息，计算出光子的能量erg。另外本文还修改了MCNP部分模块，使其能够支持光子输运计算的并行化。

1.4 并行计算

精确计算具有复杂结构的装置的停机剂量率，需要对其进行数量庞大的网格划分。为了解决计算耗时问题，本文采用多节点和多线程计算技术。

程序设计过程中，以其中一个节点作为主节点运行网格停机剂量率计算程序，负责计算文件的发送与结果文件的接收。其余的节点负责调用FISPACT进行活化计算，并将计算输出文件发送给主节点。主节点上的主程序完成计算输入文件准备之后，启动两个线程，一个用于将计算文件平均发送给计算节点进行计算；另一个负责监听计算节点传递结果文件的请求，一旦接收到计算节点的请求，主程序启动一个新的线程接收这个计算节点发回的计算结果文件。所有计算节点一直处于监听状态，一旦监听到接收文件的请求，开始接收文件。当所有文件接收完毕，关闭通信线路，并开始调用FISPACT计算。当所有文件计算完毕，计算节点再次连接主节点，并请求发送文件。

其次，在计算节点该程序创建并启动多个子线程，每个线程函数调用FISPACT进行计算，计算完成后子线程自动结束。主线程和子线程之间通过全局变量实现通信，主线程必须等待所有子线程结束后，再进行后续工作。各个子线程之间没有信息传递，从而简化了实现步骤。采用多线程的方式提高了CPU的利用率，减少了计算时间。

2 ITER-T426实验例题测试

在意大利开展的ITER-T426停机剂量率实验使用中子产生器FNG(Frascati Neutron Generator)产生的D-T中子来照射实验装置，然后使用探测器测量不同冷却时间该装置中的剂量率，将实验结果与程序计算结果进行比较分析，以校核相关计算程序、方法及核数据库的可靠性[8-10]。

ITER-T426模型由两部分组成，中子产生器模型及立体屏蔽层模型，本文使用超级蒙卡核模拟软件系统SuperMC[11-14]建立了该装置模型，如图4所示。

图4 ITER-T426模型Fig.4 ITER-T426 model

在网格计数停机剂量率计算中，对于中子输运计算，本文使用HENDL3.0[15-16]数据库计算计数网格的中子通量分布。FISPACT程序活化计算采用活化库EAF-2007[7]。中子及衰变光子输运计算使用相同的计算模型。本文在圆柱坐标系下在模型R、θ、Z方向分别划分43×8×66计数网格。在衰变光子的输运计算中，本文使用衰变光子取样程序进行衰变光子源抽样。在光子输运计算时，使用ITER推荐的通量剂量转换因子[16]计算出模型在不同停机时间的剂量率。计算结果的可视化分析使用FDS团队自主开发的科学计算可视化分析软件SVIP[18-20]。图5是该装置的中子通量分布的剖面图；图6显示了停堆12天后的光子源分布情况，可以看出右侧的不锈钢部分活化比较严重；图7展示了在图6的衰变光子源输运之后的分布情况。

图5 圆柱剖面的中子通量场(cm-2s-1)Fig.5 The profile of the neutron flux map (cm-2s-1)

图6 停机12天后中心剖面衰变光子强度分布(cm-3s-1)Fig.6 The profile of intensity of decay gamma sources 12 days after shutdown (cm-3s-1)

图7 停机12天后y-z平面光子场分布(cm-2s-1)Fig.7 The profile of photon flux map 12 days after shutdown (cm-2s-1)

图8 计算值与实验值的比值Fig.8 Comparison between calculated and experiment dose rate in the cavity center.

图8给出了网格计数停机剂量率程序计算结果与实验测量结果随冷却时间的变化。从该图中可以看出，当冷却时间比较短时，程序计算结果比实验结果最大处约低估25%。这是因为计算使用的辐照方案是简化的，不能够完全模拟实际的实验情况，并且有些对剂量贡献较多的核素没有被考虑。而冷却时间比较长时，程序计算结果最高比实验结果高出10%，此时的误差的主要来自于网格划分不够细。本程序计算结果与实验测量结果以及其他程序计算结果趋势相吻合。

3 结束语

本文根据严格二步法的计算理论，研究了基于网格计数的停机剂量率计算方法，设计并实现了基于网格计数的停机剂量率计算程序。

为了验证方法和程序的正确性，利用国际热核聚变实验堆停堆剂量基准实验ITER-T426实验进行测试，计算结果证明了该方法和程序的正确性以及处理复杂几何装置的适用性。

致谢

本工作得到中科院核能安全技术研究所·FDS团队其他成员的大力帮助和支持，在此深表感谢！

[1] Davide Valenza, Hiromasa Iida, Romano Plenteda, et al, Proposal of Shutdown Dose Rate Estimation Method by Monte Carlo Code[J]. Fusion Engineering and Design. 2001, 55: 411-418.

[2] Y. Chen, Y. Wu, U. Fischer, The rigorous 2-step calculation of shutdown dose rates for nuclear fusion devices[J]. Nuclear technology. 2003, 26: 763-764.

[3] 吴宜灿，李静惊，李莹，等. 大型集成多功能中子学计算与分析系统VisualBUS的研究与发展[J]. 核科学与工程，2007, 27(4): 365-368.

[4] L. L. Wu, Q. Yang, J. Zou, et al. Shutdown dose rate calculation code system and its application to EAST tokamak. Fusion Engineering and Design, 2012, 87: 1315-1318.

[5] 吴亮亮，应栋川，邱岳峰，等. 三维停堆剂量率计算程序研发及其在EAST上的应用[J]. 核科学与工程，2011, 31(1): 80-85.

[6] MCNP Version 5, Diagnostic applications group(X-5)[J]. Los Alamos National Laboratory.5,2003:334-567.

[7] R.A. Forrest, FISPACT2007: User Manual[M]. UKAEA Fusion. 534, 2007:12-22.

[8] P. Batistoni, M. Angelone, L. Petrizzi, et al, Experi-mental Validation of Shut-Down Dose Rates, Final Report, June 2001.

[9] P. Batistoni, M. Angelone, L. Petrizzi, et al, “Bench-mark Experiment for the validation of shut down activa-tion and dose calculation in a fusion device”, Journal of Nuclear Science and Technology, Sup. 2, p. 974-977 (August 2001), ND2001.

[10] P. Batistoni, L. Petrizzi, Task T426- Neutronics Experiments, Experimental Validation of Shut Down Dose Rates, EFF-Doc-726, March 2000.

[11] 吴宜灿，李莹，卢磊，等. 蒙特卡罗粒子输运计算自动建模程序系统的研究与发展[J]. 核科学与工程，2006, 26(1):20-27.

[12] 曾勤，卢磊，李莹，等. 蒙特卡罗粒子输运计算自动建模程序MCAM在ITER核分析建模中的应用[J]. 原子核物理评论，2006, 23(2):138-141.

[13] Y. Wu, FDS Team. CAD-based Interface Programs for Fusion Neutron Transport Simulation. Fusion Engineering and Design, 2009, 84: 1987-1992.

[14] 吴宜灿，胡丽琴，龙鹏程，等. 先进核能系统设计分析软件与数据库研发进展[J]. 核科学与工程，2010, 30(1): 55-64.

[15] J. Zou, Z. He, Q. Zeng, et al. Development and testing of multigroup library with correction of self-shielding effects in fusion-fission hybrid reactor. Fusion Engineering and Design, 2010,85: 1587-1590.

[16] C.Cheng, J. Zou, D. Xu, et al. Development and benchmark of high energy continuous-energy neutron cross section library HENDL-ADS/MC. Nuclear Science and Engineering, 2012, 32(4): 74-79.

[17] M.J. Loughlin. Recommendations on Computation of Dose from Flux Estimates. ICRP 74, 2008.

[18] Y. T. Luo, P. C. Long, G. Y. Wu, et al. SVIP-N 1.0: An integrated visualization platform for neutronics analysis. Fusion Engineering and Design, 2010, 85(7-9): 1527-1530.

[19] P. C. Long, Q. Zeng, T. He, et al. Development of a geometry-coupled visual analysis system for MCNP. Progress in nuclear science and technology, 2010, 2: 280-283.

[20] 龙鹏程, 罗月童, 邹俊, 等. 基于可编程图形处理器的可视化技术在中子学分析中的应用研究 [J]. 核电子学与探测技术, 2010, 30(8): 1042-1045.

Rigorous Two-step Shutdown Dose rate Calculation Method Based on Mesh Tally

ZHENG Jian, LI Bin, ZOU Jun, YANG Qi

(Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei, Anhui, 230031, China)

Based on the computation theory of rigorous two-step, a shutdown dose rate calculation code utilizing mesh tally technique was developed. This code supported cylindrical mesh. It integrated a source subroutine which could provide sources sampling from repeated structures geometry. Additionally, based on the platform of multi-core blade server, multi-node and multi-thread technique were adopted to conduct the mesh material and inventory calculations. The verification test was conducted on ITER-T426 experimental benchmark. The results were in good agreement with the experimental values, which proved the validity and reliability of this method and code.

Mesh tally; Shutdown dose rates; Rigorous two-step method; Automatic coupling

2016-11-21

国家ITER 973计划(2011GB11306)；中科院战略性先导科技专项(XDA03040000)；中国科学院知识创新工程重要方向项目(095CF2R211、KJCX2-YW-N35)

郑 剑(1988—)，男，湖北黄冈，硕士研究生，现主要从事中子学计算研究工作

杨 琪:qi.yang@fds.org.cn

TL63

A

0258-0918(2016)06-0784-06