1. 中国辐射防护研究院, 保健物理所剂量学室, 太原 030006;
2. 清华大学 工程物理系, 粒子技术与辐射成像教育部重点实验室, 北京 100084
收稿日期：2019-03-13
基金项目：国家自然科学基金资助项目（11305146，11675151）
作者简介：何良(1995-), 男, 博士研究生
通信作者：李君利, 教授, E-mail:lijunli@tsinghua.edu.cn
摘要：三维辐射场数据是进行剂量评估与优化的基础，为了更好地解决辐射场数据的快速计算与三维可视化显示的难题，该文基于辐射场梯度变化提出了一种变权重的非均匀网格划分方法，可根据辐射场的梯度变化划分出疏密程度不同的三维计算网格，兼顾辐射场数据的计算精度与计算效率，并使用点核积分程序对该网格划分方法进行了计算验证。结果表明，在计算精度几乎相同时，变权重的非均匀网格划分方法能够有效地缩短辐射场的计算时间，较大程度地提高点核积分方法的计算效率，同时该网格划分方法也表现出更好的辐射场显示效果。该研究可为点核积分方法在核设施三维辐射场可视化计算与剂量评估方面的应用提供技术支持。
关键词：三维辐射场点核积分网格划分剂量可视化
Mesh generation method with variable weights based on gradient changes in a radiation field
HE Liang1,2, LI Hua1, ZHAO Yuan1, LIU Liye1, CAO Qinjian1, LI Junli2
1.Department of Health Physics, Radiation Dosimetry Laboratory, China Institute for Radiation Protection, Taiyuan 030006, China;
2.Department of Engineering Physics, Tsinghua University, Key Laboratory of Particle & Radiation Imaging, Ministry of Education, Beijing 100084, China
Abstract: 3-D radiation field data is used for dose evaluations and design optimization. This paper presents a mesh generation method that uses variable weights based on gradient changes in a radiation field to improve the calculational speed and 3-D visualizations of the radiation field. The 3-D mesh has different densities that are arranged based on the radiation field gradients to balance the computational accuracy and efficiency. The point kernel integral method was used to verify the method. The results show that this non-uniform mesh generation method with variable weights effectively shortens the computing time for the point kernel integral method results. This method also provides better visualization of the radiation field. Thus, this study reduces the calculational time for the point kernel integral problem and gives better results for visualization and dose assessments of 3-D radiation fields in nuclear facilities.
Key words: 3-D radiation fieldpoint-kernel integrationmesh generationdose visualization
据国际原子能机构(IAEA)统计，截至2017年7月，全世界有447座核反应堆运行，60座在建。在高预期情况下，预计2050年核电装机量将比2016年增长123%，而中亚和东亚地区最高增长3.5倍[1]。中国作为核电大国，核电装机容量也在不断提升。国务院发布的《能源发展“十三五”规划》中显示，预计到2020年，中国核电运行容量有望达到5 800万kW，在建容量可达3 000万kW以上[2]。
核电的快速发展也给辐射防护领域带来了新的挑战。目前在核设施的运行防护中最重要的原则之一就是最优化原则(ALARA)，即考虑经济及社会等因素与安全防护水平的关系，合理分析并进行辐射防护最优化决策，尽可能地降低辐射剂量[3]。ALARA中降低剂量有多种具体措施，例如工作前模拟训练，优化工作人员的工作位置和方向[4]，以及利用模拟计算手段进行作业剂量评估与方案优化等。为此，美国、法国、日本、比利时等多个核电大国都开展了大量研究，并开发了相应软件投入实际应用[5-8]。如比利时核能研究中心(SCK·CEN)开发的VISIPLAN软件，以及法国原子能委员会(CEA)开发的NARVEOS系统等。国内在核设施三维辐射场模拟计算方面开展的研究也较多，如张永领等[9-10]基于DELMIA和VIRTOOLS平台开发的反应堆退役三维仿真原型系统，提出了三维辐射场计算和可视化技术方案；唐绍华等[11]开发了三维剂量场评价系统(RPOS)，用于核设施设计、退役及防护优化等。作业剂量模拟与优化系统的核心是核设施所在空间的三维辐射场数据，上述软件对于三维辐射场数据的计算大都使用点核积分方法作为核心算法，具有较高的计算效率，然而随着ALARA对三维辐射场的计算精度及显示效果要求的提高，辐射场计算所需三维点阵数量越来越庞大，这就对模拟计算软件整体的计算效率提出了新的挑战。如何在保证计算精度的情况下，快速进行三维辐射场的计算成为了亟待解决的问题。
三维辐射场数据的获取实际上是对网格计算点数据的计算，合理的三维网格划分会使得计算精度与计算效率之间存在较好的平衡。本文基于辐射场梯度变化提出了一种变权重的非均匀网格划分方法，并利用点核积分计算对其进行了分析，验证了该方法在一定范围内的适用性，旨在为辐射场可视化计算及模拟仿真领域提供技术支持。
1 材料和方法三维辐射场是包含坐标和剂量信息的一系列三维点阵，通过网格划分实现。在实际应用中，辐射场网格的精细程度对剂量计算及可视化显示影响很大。一方面从辐射场可视化的角度，原始辐射场网格往往间距较大，不能满足可视化显示效果的要求，需要基于原始的网格数据通过插值计算得到更为精细的纹理矩阵；另一方面从剂量计算角度，需要根据人员的路径进行剂量计算，同样需要对原始网格进行细化。辐射场网格的进一步细化如图 1所示。
 |
| 图 1 辐射场网格的进一步细化 |
| 图选项 |
辐射场网格划分越精细，所需计算的点个数越多，计算时间也就越长。当体源或面源的离散数目一定时，计算点的个数与计算时间成正比，并且当网格数量达到一定程度后，计算精度会趋于饱和。如果对整个计算区域均采用尺寸较小的网格划分，虽然提高了计算精度，但同时会造成计算时间的大幅增加，并且对于辐射场变化不大的区域，细小的网格划分对计算精度提升并不高。综合考虑计算精度与时间，本文在计算辐射场时根据辐射场变化梯度不同采用非均匀网格划分，即在靠近辐射源的区域，由于辐射场变化梯度较大可划分较细的网格，在辐射场变化缓慢的区域则采用较粗的网格。目的是在保证基于辐射场数据的作业剂量计算结果精度的前提下，尽可能地提高辐射场计算的效率。
不同形状的辐射源在其附近所产生的辐射场的梯度变化也有所不同。以线源为例，如图 2所示，取线源中垂线上的计算点，并且不考虑屏蔽情况，对其辐射场变化梯度进行推导。
 |
| 图 2 长度为L的线状放射源 |
| 图选项 |
设γ线源的长度为L，单位长度的活度为A。计算点与线源距离为r，x代表线源上点的位置。则线源中一小段dx可以看作为点源，其活度为Adx。根据点源计算公式，可得到该小段源在计算点处的剂量率，并对其从-L/2至L/2进行积分，得到整个线源在计算点处的剂量率D为
| $D = \int_{ - L/2}^{L/2} {\frac{{\mathit{F} \cdot \mathit{A}}}{{({r^2} + {x^2})}}} dx = \frac{{2AF}}{r}{\rm{t}}{{\rm{g}}^{ - 1}}\frac{L}{{2r}}.$ | (1) |
| $\frac{{\partial D}}{{\partial r}} = - \frac{{2AF}}{r}{\rm{t}}{{\rm{g}}^{ - 1}}\frac{L}{{2r}} - \frac{{AFL}}{{{r^3} + {L^2}r/4}}.$ | (2) |
 |
| 图 3 不同形状放射源其附近辐射场的变化 |
| 图选项 |
针对上述问题，本文采用变权重的非均匀网格划分的方式来确定计算点，具体分为4个步骤，如图 4所示。
 |
| 图 4 变权重非均匀网格划分的流程图 |
| 图选项 |
1) 首先根据实际需求自定义设置网格参数进行均匀划分，考虑到辐射场数据是用于作业人员受照剂量计算的目的，最大网格间距参考正常人体肩宽的距离，一般不宜超过0.5 m；
2) 针对不同类型的源项，需选取不同的权重倍数对原始网格进行细化处理。权重倍数的选择与源的类型有关，既要考虑源附近辐射场的梯度变化，还要结合实际需要，综合考虑。一般而言，对于点源、线源附近区域采用的权重倍数为3~5，即在其细化区域采用的网格间距为原始的1/3~1/5；对于球源、圆柱体源等其他类型的源项附近区域，采用的权重倍数为2~3，即在其细化区域采用的网格间距为原始的1/2~1/3。
3) 在源位置所在的网格及相邻网格中，包含前后、左右、上下3个方向，分别对应X、Y、Z轴。采用步骤2中选取的权重倍数，对原始网格进行细化。这里所指的相邻网格是指体源或面源外表面所在网格的相邻网格，在实际计算中，细化区域一般向外延伸1个原始网格就能满足大部分情况，在少数情况下，可以向外延伸2~3个原始网格。以点源为例，其非均匀网格划分的过程如图 5所示，这里仅给出了二维的网格细化过程。
 |
| 图 5 点源附近非均匀网格划分的过程 |
| 图选项 |
4) 如果计算或测量过程中存在多个源项，重复步骤2和3，直至对于所有的源项附近区域，均进行了网格细化处理，形成最终的三维非均匀计算网格。
2 结果与讨论为了验证本文提出的非均匀网格划分方法的实用性，以某实验室房间内平面辐射场的计算为例，对传统的均匀网格划分和本文的非均匀网格划分方法进行对比分析，该房间的俯视图如图 6所示。
 |
| 图 6 某房间内的俯视图 |
| 图选项 |
该实验室房间长570 cm、宽870 cm、高450 cm，房间内有一管道源项，管道外半径10 cm、内半径8 cm，管壁材料为铁，其中心轴线与Y轴平行，位于x=490 cm，z=370 cm处，其内表面均匀分布一层单能的γ面源，总活度为8×108 Bq，γ射线的能量为1.173 MeV。为了便于比较，在计算辐射场时，只选取高度为360 cm处的平面网格。
在本算例中，γ剂量率的计算采用的是点核积分方法。点核积分方法是一种常用的剂量模拟计算方法，它在处理复杂几何空间的辐射屏蔽问题效果很好。相比于Monte Carlo方法，点核积分原理简单，计算快速，适用于进行γ射线的屏蔽计算[12]。本文计算所用的点核积分程序由中国辐射防护研究院自主研发，采用了较新的累积因子数据[13-14]，具有较高的计算精度与效率。
本算例在辐射场网格划分时，采用4种不同的网格进行，分别如下所示：
1) 采用13×20的均匀网格进行划分；
2) 采用40×60的均匀网格进行划分；
3) 首先采用13×20的均匀网格进行划分，然后在对源区所在网格及其相邻网格均采用权重倍数为2进行细化；
4) 首先采用13×20的均匀网格进行划分，然后在对源区所在网格及其相邻网格均采用权重倍数为3进行细化。
在同一台计算机上，对上述4种不同网格划分下的平面辐射场进行计算，计算时管道面源的离散数目均不变，所用的计算时间如表 1所示。
表 1 不同网格划分下的辐射场计算所用的时间
| 网格划分 | 13×20 | 40×60 | 局部2倍细化 | 局部3倍细化 |
| 计算时间/s | 6 | 61 | 12 | 20 |
表选项
为了对比这4种不同网格划分下的剂量计算的准确性，在不同的辐射场数据下，分别采用线性插值计算一组位置处的剂量率，进行对比。为了方便，做如下标记：
基于13×20划分下的网格数据插值计算得到的某位置处的剂量率记为D1；40×60划分下剂量率的数值记为D2；同理局部3倍加密下剂量率的数值记为D3；局部2倍加密下剂量率的数值记为D4；同时点核积分程序在该位置处直接计算的剂量率数值记为D。
在上述不同网格划分下，得到的不同位置处剂量率结果的对比如表 2所示。可以看出，对于距离源较远的点，4种网格划分下线性插值得到的辐射场数据均比较接近实际值；对于距离源较近的点，显然网格划分得越细，得到的结果越接近实际值，网格较粗时插值计算会产生较大的误差。同时在本算例中，局部3倍细化网格下的辐射场计算精度与40×60网格下的计算精度几何相同，但前者所用的计算时间仅是后者所用时间的1/3。
表 2 不同网格划分下辐射场数据的对比
| 位置坐标/cm | 剂量率/(μGy·h-1) | |||||||
| X | Y | Z | D1 | D2 | D3 | D4 | D | |
| 50 | 50 | 360 | 2.20 | 2.20 | 2.20 | 2.21 | 2.20 | |
| 100 | 100 | 360 | 2.85 | 2.84 | 2.85 | 2.85 | 2.84 | |
| 200 | 200 | 360 | 4.99 | 4.97 | 4.99 | 4.99 | 4.97 | |
| 300 | 300 | 360 | 9.35 | 9.22 | 9.35 | 9.35 | 9.20 | |
| 400 | 400 | 360 | 22.37 | 20.89 | 22.37 | 22.37 | 20.80 | |
| 450 | 450 | 360 | 67.32 | 46.82 | 47.74 | 46.47 | 46.00 | |
| 460 | 460 | 360 | 83.99 | 63.44 | 61.44 | 70.09 | 60.30 | |
| 470 | 470 | 360 | 100.66 | 92.15 | 93.15 | 96.03 | 86.40 | |
表选项
此外，在辐射场显示方面，本文中非均匀网格划分方法也大大改善了源附近区域的辐射场显示效果。如图 7所示，将算例中的平面辐射场进行可视化显示，对比两种不同网格划分下的显示效果，显然图 7b比图 7a更加对称且光滑，与理论相一致，效果更好。
 |
| 图 7 (网络版彩图)不同网格划分下辐射场显示效果的对比 |
| 图选项 |
3 结论三维网格点阵会对辐射场可视化计算的精度及效率产生非常大的影响。本文基于辐射场梯度变化提出了一种变权重非均匀网格计算点的划分方法，即在传统均匀网格划分的基础上，采用权重倍数进行局部网格细化，对辐射场梯度变化大的区域采用较小的网格间距，而对辐射场数值小或梯度变化小的区域采用较大的网格间距。同时以某实验室房间内平面辐射场的计算为例，使用点核积分方法对该网格划分方法进行了验证。
结果表明，相比传统的均匀网格划分，在计算精度几乎相同的条件下，本文方法能够有效地缩短辐射场计算的时间，较大程度地提高点核积分方法计算三维辐射场的效率，同时大大改善了放射源附近区域的显示效果，可为核设施现场的辐射防护最优化提供先进的技术支持手段。
参考文献
| [1] | IAEA. International status and prospects for nuclear power 2017[R/OL]. (2017-07-28)[2019-01-10]. https://www-legacy.iaea.org/About/Policy/GC/GC61/GC61InfDocuments/English/gc61inf-8_en.pdf |
| [2] | 王丹. 解读"十三五"发展规划, 助力核电发展[J]. 中国核电, 2017(2): 157-160. WANG D. Interpretation of the 13th five-year development plan in support of nuclear power development[J]. China Nuclear Power, 2017(2): 157-160. (in Chinese) |
| [3] | 刘华. 辐射防护最优化方法及其应用[J]. 核安全, 2007(2): 1-6. LIU H. Optimization methods of radiation protection and its application[J]. Nuclear Safety, 2007(2): 1-6. DOI:10.3969/j.issn.1672-5360.2007.02.001 (in Chinese) |
| [4] | 潘志强. 辐射安全手册[M]. 北京: 科学出版社, 2011. PAN Z Q. Radiation safety manual[M]. Beijing: Science Press, 2011. (in Chinese) |
| [5] | SAUNDERS P, RAHON T, QUINN D, et al. Demonstration of advanced 3-D ALARA planning prototypes for dose reduction[R]. California: Electric Power Research Institute (EPRI), 2012. |
| [6] | MóL A C A, PEREIRA C M N A, FREITAS V G G, et al. Radiation dose rate map interpolation in nuclear plants using neural networks and virtual reality techniques[J]. Annals of Nuclear Energy, 2011, 38(2-3): 705-712. DOI:10.1016/j.anucene.2010.08.008 |
| [7] | YUKIHARU O, MITSUKO F, KIYOTAKA S, et al. A System for the calculation and visualisation of radiation field for maintenance support in nuclear power plants[J]. Radiation Protection Dosimetry, 2005, 116(1-4): 592-596. DOI:10.1093/rpd/nci014 |
| [8] | VERMEERSCH F. ALARA pre-job studies using the VISIPLAN 3D ALARA planning tool[J]. Radiation Protection Dosimetry, 2005, 115(1-4): 294-297. DOI:10.1093/rpd/nci115 |
| [9] | 戴波, 张永领, 周斌, 等. 反应堆退役仿真技术研究方案[J]. 核动力工程, 2013, 34(3): 168-171. DAI B, ZHANG Y L, ZHOU B, et al. Study scheme of reactor decommissioning simulation technology[J]. Nuclear Power Engineering, 2013, 34(3): 168-171. DOI:10.3969/j.issn.0258-0926.2013.03.040 (in Chinese) |
| [10] | 张永领, 胡一非, 刘猛, 等. 反应堆退役三维辐射场实时计算及可视化[J]. 辐射防护, 2018(1): 19-25. ZHANG Y L, HU Y F, LIU M, et al. Real-time computation and visualization of 3-D radiation field for nuclear reactor decommissioning scene[J]. Radiation Protection, 2018(1): 19-25. DOI:10.3969/j.issn.1004-6356.2018.01.005 (in Chinese) |
| [11] | 唐邵华, 吕炜枫, 刘杰, 等. 核电站三维剂量场评价系统的开发及应用[J]. 辐射防护, 2017, 37(5): 347-354. TANG S H, LV W F, LIU J, et al. Development and application of 3-D dose rate field evaluation system[J]. Radiation Protection, 2017, 37(5): 347-354. (in Chinese) |
| [12] | 李春槐, 张立吾. 几何空间装配法在点核积分计算中的应用[J]. 核动力工程, 1988(5): 42-52. LI C K, ZHANG L W. Application of geometry space configuration method in point-kernel integral calculation[J]. Nuclear Power Engineering, 1988(5): 42-52. (in Chinese) |
| [13] | 李华, 赵原, 刘立业, 等. 基于MCNP对γ射线吸收剂量累积因子的计算与研究[J]. 辐射防护, 2017, 37(3): 161-168. LI H, ZHAO Y, LIU L Y, et al. Research on gamma ray buildup factor for energy absorption based on MCNP[J]. Radiation Protection, 2017, 37(3): 161-168. (in Chinese) |
| [14] | 李华, 赵原, 刘立业, 等. 介质尺寸对水中γ射线吸收剂量累积因子的影响[J]. 清华大学学报(自然科学版), 2017, 57(5): 525-529. LI H, ZHAO Y, LIU L Y, et al. Effect of medium size on the γ-ray buildup factor for energy absorption in water[J]. Journal of Tsinghua University(Science and Technology), 2017, 57(5): 525-529. (in Chinese) |
