据国际原子能机构(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所示。

辐射场网格划分越精细，所需计算的点个数越多，计算时间也就越长。当体源或面源的离散数目一定时，计算点的个数与计算时间成正比，并且当网格数量达到一定程度后，计算精度会趋于饱和。如果对整个计算区域均采用尺寸较小的网格划分，虽然提高了计算精度，但同时会造成计算时间的大幅增加，并且对于辐射场变化不大的区域，细小的网格划分对计算精度提升并不高。综合考虑计算精度与时间，本文在计算辐射场时根据辐射场变化梯度不同采用非均匀网格划分，即在靠近辐射源的区域，由于辐射场变化梯度较大可划分较细的网格，在辐射场变化缓慢的区域则采用较粗的网格。目的是在保证基于辐射场数据的作业剂量计算结果精度的前提下，尽可能地提高辐射场计算的效率。

不同形状的辐射源在其附近所产生的辐射场的梯度变化也有所不同。以线源为例，如图 2所示，取线源中垂线上的计算点，并且不考虑屏蔽情况，对其辐射场变化梯度进行推导。

设γ线源的长度为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)

其中F为转换系数，只与γ光子的能量有关。为了得到该剂量率随r的梯度变化，对式(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)

其中A、L、F在放射源的信息确定后不会再发生改变，即为常量。辐射场的梯度变化会随着r的增大，呈现约立方衰减的趋势，并逐渐趋向于0。其他形状的放射源产生的辐射场变化具有类似的趋势，如图 3所示。图中给出了点源、线源、圆柱体源，以及球型体源产生的辐射场变化趋势，图中纵坐标为归一化的γ剂量率，横坐标采用对数刻度，为计算点到源表面的距离，线源和圆柱源计算点取自过其中心的与轴垂直方向上的点。可以看出，点源附近的辐射场变化程度最大，遵循平方反比的关系；线源长度和球源的半径越短，其附近的辐射场变化趋势越接近点源的情况；圆柱体源越细，附近的辐射场变化趋势越接近线源的情况。

针对上述问题，本文采用变权重的非均匀网格划分的方式来确定计算点，具体分为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所示，这里仅给出了二维的网格细化过程。

4) 如果计算或测量过程中存在多个源项，重复步骤2和3，直至对于所有的源项附近区域，均进行了网格细化处理，形成最终的三维非均匀计算网格。

2 结果与讨论

为了验证本文提出的非均匀网格划分方法的实用性，以某实验室房间内平面辐射场的计算为例，对传统的均匀网格划分和本文的非均匀网格划分方法进行对比分析，该房间的俯视图如图 6所示。

该实验室房间长570 cm、宽870 cm、高450 cm，房间内有一管道源项，管道外半径10 cm、内半径8 cm，管壁材料为铁，其中心轴线与Y轴平行，位于x=490 cm，z=370 cm处，其内表面均匀分布一层单能的γ面源，总活度为8×10⁸ 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所示。

为了对比这4种不同网格划分下的剂量计算的准确性，在不同的辐射场数据下，分别采用线性插值计算一组位置处的剂量率，进行对比。为了方便，做如下标记：

基于13×20划分下的网格数据插值计算得到的某位置处的剂量率记为D₁；40×60划分下剂量率的数值记为D₂；同理局部3倍加密下剂量率的数值记为D₃；局部2倍加密下剂量率的数值记为D₄；同时点核积分程序在该位置处直接计算的剂量率数值记为D。

在上述不同网格划分下，得到的不同位置处剂量率结果的对比如表 2所示。可以看出，对于距离源较远的点，4种网格划分下线性插值得到的辐射场数据均比较接近实际值；对于距离源较近的点，显然网格划分得越细，得到的结果越接近实际值，网格较粗时插值计算会产生较大的误差。同时在本算例中，局部3倍细化网格下的辐射场计算精度与40×60网格下的计算精度几何相同，但前者所用的计算时间仅是后者所用时间的1/3。

此外，在辐射场显示方面，本文中非均匀网格划分方法也大大改善了源附近区域的辐射场显示效果。如图 7所示，将算例中的平面辐射场进行可视化显示，对比两种不同网格划分下的显示效果，显然图 7b比图 7a更加对称且光滑，与理论相一致，效果更好。

3 结论

三维网格点阵会对辐射场可视化计算的精度及效率产生非常大的影响。本文基于辐射场梯度变化提出了一种变权重非均匀网格计算点的划分方法，即在传统均匀网格划分的基础上，采用权重倍数进行局部网格细化，对辐射场梯度变化大的区域采用较小的网格间距，而对辐射场数值小或梯度变化小的区域采用较大的网格间距。同时以某实验室房间内平面辐射场的计算为例，使用点核积分方法对该网格划分方法进行了验证。

结果表明，相比传统的均匀网格划分，在计算精度几乎相同的条件下，本文方法能够有效地缩短辐射场计算的时间，较大程度地提高点核积分方法计算三维辐射场的效率，同时大大改善了放射源附近区域的显示效果，可为核设施现场的辐射防护最优化提供先进的技术支持手段。