三维停堆剂量率计算程序研发及其在EAST上的应用

吴亮亮,应栋川,邱岳峰,王国忠,张延云,宋 婧何 桃,,周少恒,熊 健,,李 佳,

(1.中国科学技术大学核科学技术学院,安徽 合肥 230027;2.中国科学院等离子体物理研究所,安徽合肥 230031;3.合肥工业大学数学学院,安徽 合肥 230069)

停堆剂量是由核装置运行期间被中子活化的材料所发射的衰变光子辐射引起的。在核装置在运行停止或运行间歇期间,工作人员需要靠近或进入装置内部进行实验测量或维修检测工作,将受到放射性照射的危险。因此准确地评估核装置的停堆剂量率水平,对于核装置设计及维护人员的安全有着重要的参考意义[1]。

对于聚变装置停堆剂量率的计算,通常采用的是基于三维蒙特卡罗程序的停堆剂量率计算方法[2-3]。由于该方法仅能考虑有限的核素,且只能处理一步的衰变过程,不适用于求解具有复杂的几何结构先进反应堆(如核聚变反应堆、聚变裂变混合堆、加速器驱动次临界堆等)的停堆剂量率。本文在FDS团队自主研发的大型集成多功能中子学计算与分析系统VisualBUS[4]框架下,通过采用严格二步法[1,5](R2S：Rigorous 2-step method),耦合具有处理复杂几何能力的三维蒙特卡罗程序MCNP[6]和国际上广泛使用的活化计算程序FISPACT[7],发展了三维停堆剂量计算程序。该程序实现了中子输运、材料活化和光子剂量计算的自动耦合,可获得核装置停堆后周围空间剂量场的三维精细分布信息。本文将从计算方法、程序实现及其在先进实验超导托卡马克EAST(Experimental Advanced Superconducting Tokamak)上的初步应用三个方面来介绍该程序的发展状况。

1 停堆剂量率计算方法

核装置在运行期间,部分组件在中子的辐照之下被活化而具有放射性。装置停止运行或运行间歇期间,被活化了的材料会向周围空间发出γ光子[8]。因此,要计算停堆剂量率的分布情况,需要先求得衰变光子源在时间、空间及能量上的分布。R2S方法计算过程如图1所示。

图1 R2S方法计算流程Fig.1 Flow scheme of R2Smethod

对于合适的计算模型,中子输运计算可以得到装置运行过程中所有非空栅元中子注量率在空间、能量上的分布情况。对于几何结构复杂的三维模型(如托卡马克装置),中子输运计算通常采用适合求解复杂三维几何问题的蒙特卡罗方法[9]。本文利用国际上通用的三维蒙特卡罗输运程序MCNP进行中子输运计算以获得下一步活化计算所需的中子注量率数据。

R2S方法利用中子输运计算所得的中子注量率数据,在一定的辐照时间和冷却时间的条件下,进行材料活化计算。本文使用国际上广泛使用的活化计算程序FISPACT-2007得到光子源的描述文件,用于下一步的光子输运计算。

最后,R2S方法进行光子输运模拟,得到不同冷却时间下核装置周围空间的三维剂量场分布。计算中所使用的数据库为国际原子能机构(IAEA)最新版本聚变评价数据库FENDL-2.0[10]。

2 程序实现

基于R2S方法的三维停堆剂量计算程序采用两个接口程序将蒙特卡罗输运计算程序MCNP和活化计算程序FISPACT自动连接,实现精确求解核装置的停堆剂量率。三维停堆剂量率计算程序流程如图2所示。

其中接口程序 MF-Interface的主要功能是解析输运计算的输出文件,自动生成包含中子注量率和材料信息的活化计算输入文件;接口程序FM-Interface的主要功能是处理活化计算输出文件,将衰变光子强度及能谱分布计算结果经处理后生成光子源的描述文件,用于MCNP衰变光子输运计算。

图2 停堆剂量率计算程序流程图Fig.2 Shutdow n dose calculation code based on R2S

为提高计算效率,程序允许用户提供输运计算输出文件,避免重复进行MCNP程序计算。同时,根据用户的不同需求,计算中可以采用多种中子能群结构(包括69群的W IMSD、175群的VITAMIN-J和 315群的 TRIPOLI等能群结构)[7]。

在实际计算中,组成栅元材料的某些微量元素在输运计算中可以被忽略,但可能对活化的结果造成巨大的影响,考虑这些微量元素对提高活化计算的精确度具有重要意义。为此,程序提供了用户自定义材料的功能,用户可以更改活化计算的材料,加入输运中被忽略的微量元素,提高了活化计算的精确程度。

由于输运计算程序MCNP对源卡个数的限制,无法进行复杂的核装置模型的计算。本文通过使用接口程序FM-Interface对栅元与源卡的对应关系进行优化,使得源卡的个数满足MCNP程序的限值。同时,程序提供一个方便的源描述接口,用户可以通过定义一个包围所有栅元的较大取样空间来完成对衰变光子源的随机取样,从而使MCNP衰变光子输运计算成为可能。

3 在EAST上的初步应用

3.1 计算模型

EAST是中科院等离子体物理研究所自主设计建造的世界上首个全超导托卡马克装置。主机部分高11m,直径8m,重400 t,由超高真空室、纵场线圈、极向场线圈、内外冷屏、外真空杜瓦、支撑系统等六大子系统组成。

本文采用未经简化EAST装置的45°CAD模型作为分析对象,通过FDS团队自主研发的蒙特卡罗计算自动建模软件MCAM[12-15]将该CAD模型转换为三维中子学MCNP计算模型(如图3所示)。该计算模型包含了精细的组件结构,使得剂量计算结果的更加精确。EAST中子学模型的材料组成如表1所示。

图3 带大厅的EAST装置45°模型Fig.3 The 45°modelof EASTw ith a hall

表1 EAST中子学模型的材料组成[15]Table1 Materia l composition of EAST neutronicsmodel

本文计算分析了EAST装置在D-D等离子体反应模式下的停机剂量信息。假设中子源均匀分布在等离子体区,能量为2.45 MeV,典型D-D反应的中子产额为1×1015n/s,使用单脉冲放电时间1 000 s作为辐照参数[16]。为了较为全面地考查停机后不同时间段的辐射剂量情况的变化趋势,本文在计算方案中选取了6个冷却时间,分别为10m in,100min,400min,1 000 m in和10 000 min。

3.2 计算分析

使用FDS团队自主研发的科学计算可视化分析软件SV IP[16]对三维停机剂量程序计算所得结果进行可视化,得到EAST装置停机后周围三维空间剂量场的分布。图4分别显示了以装置为中心的水平横向截面和竖直纵向截面上剂量场的分布情况。

图4 SV IP中EAST装置三维剂量场分布Fig.4 Th ree-dimensional dose distribution of EAST

1)剂量率在径向上的变化

分析计算不同冷却时间下托克马克窗口高度处剂量率,得到剂量率随径向变化情况的分布曲线,如图5所示。由于在靠近等离子体的区域,受到的中子辐照较强,因而辐射剂量率相对较高。从图中可以看出,在靠近托克马克表面(450 cm)处辐射剂量值最高。而随着半径方向距离的增加,辐射剂量逐渐减小。不同的冷却时间的情况下,剂量率沿径向的分布近似遵循指数衰减规律。当冷却时间超过1 000m in,装置周围的剂量水平低于 10μSv/h,满足EAST装置的辐射防护标准[15]。

图5 剂量率径向分布曲线Fig.5 The curve of dose rates along the radius

2)剂量率在极向上的变化

分析计算外真空杜附近极向上的数据对剂量率的分布情况,得到变化曲线如图6所示。可以看出,不同冷却时间的剂量率沿极向方向的分布规律大致相同,在高度约为550 cm处剂量率达到最大值,离该位置越远则剂量率越小。因为此位置最靠近堆芯等离子体区域,结构材料受到的中子辐照主要集中在该高度附近,活化光子源也主要集中在该区域。同时可以看出,在冷却时间超过1 000 m in后,装置周围的剂量水平低于10μSv/h,满足EAST装置的辐射防护标准。

3)剂量率随冷却时间的变化

图6 剂量率极向分布曲线Fig.6 The curve of dose rates along the axis

为了研究EAST装置的停机剂量率随冷却时间的变化情况以便确定剂量率水平在多长的冷却时间后可以达到对人体安全的水平,本文选取了不同冷却时间下大厅中剂量率最高处(R=450 cm,Z=550 cm,θ=45°)的剂量率进行了比较分析,得到剂量率随冷却时间的变化曲线,如图7所示。

图7 剂量率随冷却时间的变化曲线Fig.7 The curve o f dose rate with the coo ling time

由图可以看出停机剂量率在大约100 m in以前,衰减比较慢。因为这段时间内,停堆剂量率主要来自于装置材料SS304L中的杂质55Mn吸收中子发生(n,γ)俘获反应生成的放射性核56M n,半衰期为2.58 h。随后,51Cr(半衰期为27.8 d)开始变成主要放射性核素,它主要来自于装置材料SS304L中的杂质50Cr与中子发生(n,g)反应,从图中可以看出,在100 min至10 000 min(约一星期)衰减加快。另外,放射性核素55Fe、60Co及60Com在相应的冷却时间段对停堆剂量亦有所贡献。根据EAST装置的辐射防护标准,EAST装置正常运行时工作场所的辐射水平不得高于10μSv/h。从图中可以看到,当冷却时间大于 1 000 m in以后,EAST装置的剂量水平小于10μSv/h。因此在这种情况下,工作人员在停机以后必须等待1 000 m in之后,才能靠近装置进行操作。而一周以后,大厅内的剂量率则已经衰减到天然本底的水平(0.1μSv/h)。

4 结论

本文基于严格二步法(R2S),耦合蒙特卡罗输运计算程序MCNP和活化计算程序FISPACT,发展了三维停机剂量计算程序。该程序充分利用了MCNP程序中子、光子输运计算以及处理复杂几何结构的能力,并结合FISPACT程序材料活化及光子衰变计算的能力,实现了核装置停机后周围空间剂量场分布的精确描述和分析。

该停机剂量计算程序已初步应用于EAST装置,所得结果直观显示了 EAST装置停机后,不同冷却时间和空间下剂量场的三维空间分布情况,这些数据可以为EAST辐射防护工作提供了重要参考。同时,停机剂量计算程序在EAST装置上的成功应用也为其在其他核聚变装置上的应用提供了参考。在下一步的工作中,需要进行各种单脉冲放电和多脉冲放电条件下停机剂量率的分析计算,以获得EAST装置停机后辐射剂量场在三维空间的真实分布。

5 致谢

本文工作是在中科院知识创新工程重要方向项目(编号k JCX2-YW-N37与KJCX 2-YWN35)和中科院信息化专项项目(编号INFO-115-C01-SDB2-010)的资助下进行的。

[1] 陈义学,吴宜灿,Fischer U.核聚变装置停机剂量率分析计算的严格两步(R2S)法[J].核技术,2003,26(10)：763-764.

[2] Iida H,Valenza D,Plentada R,et al.JNucl Sci Tech,2000(Supp l1)：235-242.

[3] Valenza D,Iida H,Plenteda R,et al.Fus Eng Design,2001,55：411-418.

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

[5] Davis A,Pampinb R.Benchmarking theMCR2S system for high-resolu tion activation dose analysis in ITER[J].Fusion Engineering and Design,FUSION-5010,2009：1-2.

[6] Briesmeister JF,et al.MCNP-A General Mon te Carlo N-Particle T ransport Code[M].Version 4C,Los A lamos National Laboratory,Report LA-13709-M,Ap ril 2000.

[7] Forrest R A.FISPACT 2007：User Manual[R].UKAEA Fusion,Report UKAEA FUS 534,March 2007.

[8] 柴竹新,吴宜灿,刘伯学.核聚变装置 EAST高可靠性辐射防护控制系统[J].核电子学与控测技术,2005,25(1)：28-29.

[9] 谢仲生,邓力.中子输运理论数值计算方法[M].陕西：西北工业大学出版社,2005.

[10]Pashchenko A B,Wienke H,Kopecky J,et al.FENDL/A-2.0 Neutron Activation Cross Section Data Lib rary for Fusion Application[R].IAEA Vienna,Report IAEA-NDS-173,Rev.1,October 1998.

[11]W u Y C.CAD-based interface p rog ram s for fusion neutron transport simu lation[J].Fusion Engineering and Design,2009,84(7-11)：1987-1992.

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

[13]曾勤,卢磊,李莹,等.蒙特卡罗粒子输运计算自动建模程序MCAM在ITER核分析建模中的应用[J].核科学与工程,2006,23(2)：138-141.

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

[15]陈义学.托卡马克装置三维辐射场计算方法的发展及其在HT-7U环境影响评价中的应用[D].中国科学院等离子体物理研究所,2002.

[16]罗月童,龙鹏程,薛晔,等,面向中子学分析的集成可视化平台SV IP的发展研究[J].核科学与工程,2007,27(4)：374-378.