CN107194103B

CN107194103B - 一种基于3dsMax的核设施模型辐射场剂量仿真方法

Info

Publication number: CN107194103B
Application number: CN201710421231.3A
Authority: CN
Inventors: 刘永阔; 杨立群; 彭敏俊; 李梦堃; 白璐
Original assignee: Harbin Engineering University
Current assignee: Harbin Engineering University
Priority date: 2017-06-07
Filing date: 2017-06-07
Publication date: 2020-11-10
Anticipated expiration: 2037-06-07
Also published as: CN107194103A

Abstract

本发明提供的是一种基于3dsMax的核设施模型辐射场剂量仿真方法。用3dsMax软件根据确定的核设施参数构建模型，并将文件保存为3DS格式；导入3DS核设施模型文件，获得模型参数；用八叉树法将3dsMax核设施模型分解为体素；将确定的体素参数与材质信息写成输入卡；将输入卡导入到点核积分程序内；计算累积因子；计算伽马射线在辐射场中的平均自由程；利用SQLite数据库引擎建立通量率‑剂量率转换因子、各化学元素与材料的质量衰减因子以及单层累积因子数据库；运用布尔连接运算符对箱体进行组合运算，构造复杂的辐射场几何结构；用点核积分方法计算三维辐射场剂量。本发明能实现对有尺寸、材质与能量参数的复杂3dsMax核设施模型的辐射场剂量计算。

Description

一种基于3dsMax的核设施模型辐射场剂量仿真方法

技术领域

本发明涉及的是一种核退役仿真方法，具体地说是一种用3dsMax对核设施辐射场剂量仿真的方法。

背景技术

核设施维修与退役是核设施生命周期的重要环节。在制定核设施维修或退役策略时，需要优化维修、退役方案以及工作人员在退役作业过程中的路径，进而降低辐射对工作人员的伤害，这就需要准确地了解核设施三维空间辐射剂量的分布情况，进行辐射仿真。在辐射仿真过程中，为获得准确的虚拟辐射场剂量计算结果，就必须对复杂模型进行辐射计算。

目前，核退役虚拟辐射场的模拟与研究大多使用国外成熟的仿真软件来完成，其中蒙特卡罗方法与点核积分方法在屏蔽设计中有着广泛的应用。点核积分方法通过引入累积因子来考虑散射光子对辐射量的影响，将辐射场中所有源项按照几何尺寸离散为点源，并将源项能谱离散为若干离散值，然后分别计算不同能量与不同点源在各个剂量点的剂量值，最后将同一剂量点的剂量值叠加计算出探测点总的剂量值。但这些软件对复杂模型的辐射场剂量计算较为简单，无法直接获得复杂3dsMax核设施模型的精确辐射场剂量。

综上所述，开发出一种准确、可靠的计算复杂3dsMax核设施模型辐射场剂量的仿真方法对核退役仿真具有重大的实际意义。

发明内容

本发明的目的在于提供一种准确、可靠的基于3dsMax的核设施模型辐射场剂量仿真方法。

本发明的目的是这样实现的：

(1)用3dsMax软件根据确定的核设施参数构建模型，并将文件保存为3DS格式；

(2)导入3DS核设施模型文件，获得模型参数；

3DS文件由许多块(chunk)组成，每个块首先描述其信息类别。块的前两项信息分别是：块的ID和块的长度。块的ID作为块的标示，块的长度表示下一个块相对于该块起始位置的偏移字节数。3DS文件最开始出现的主块是基本块，它包含了整个文件。3D编辑程序块ED—IT3DS主要定义物体的形体数据，块中包含了一个物体描述子块EDIT_OBJECT，其下有三角形列表子块OBJ_TRIMESH。OBJ_TRIMESH主要包括了体素化程序所需的模型表面信息。

(3)用八叉树法将3dsMax核设施模型分解为体素；

本发明对模型的三维网格进行体素化获得体素模型。由于真实的大型核设施多为空心的壳结构，所以算法可以简化为对模型表面的体素化，即对组成模型表面所有三角形面的体素化。由于三角形尺寸一般都小于体素尺寸，所以算法可以简化为对所有三角形顶点与中点的体素化。

(4)将确定的体素参数与材质信息写成输入卡；

输入卡由6部分组成，包括曲面卡、基本体数据卡、材料数据卡、区域卡、放射源卡以及探测点卡。

(5)将输入卡导入到点核积分程序内；

(6)计算累积因子；

本发明采用ANSI/ANS-6.4.3数据库以及G-P拟合公式计算单层累积因子，采用经验公式计算双层或多层材料累积因子。

(7)计算伽马射线在辐射场中的平均自由程；

平均自由程计算公式为

式中，t(E)是伽马光子从点源到探测点穿过所有屏蔽材料的平均自由程；i为伽马射线穿过的空间区域编号；ρ_i为空间区域i的材料密度；μ_i(E)/ρ为在光子能量为E时，空间区域i的材料的质量减弱系数；d_i为伽马射线在区域i中的几何距离。

(8)利用SQLite数据库引擎建立通量率-剂量率转换因子、各化学元素与材料的质量衰减因子以及单层累积因子数据库；

(9)运用布尔连接运算符对箱体进行组合运算，构造复杂的辐射场几何结构；

(10)用点核积分方法计算三维辐射场剂量。

点核积分方法的基本思想是将放射源离散为若干点核，然后计算每个点核在探测点的剂量值，最后将所有点核的剂量值求和得到探测点总的剂量率。本发明中假设核反应堆以及其他放射源能够由一系列离散的各向同性点源组成，并且可以通过对各离散点核的贡献求和来获取放射源在探测点总的辐射效应。在点核积分方法的几何模型中，各点核在探测点的剂量值为

式中，r_p与r_d分别为点核与探测点的位置；E为光子能量；C(E)为伽马光子辐射效应转换因子；S(E,r_d)是点核源项强度；B(E,t)是累积因子；t(E)是伽马光子从点源到探测点穿过所有屏蔽材料的平均自由程，计算公式为

式中，i为伽马射线穿过的空间区域编号；ρ_i为空间区域i的材料密度；μ_i(E)/ρ为在光子能量为E时，空间区域i的材料的质量减弱系数；d_i为伽马射线在区域i中的几何距离。

将探测点的剂量值在整个源项体积与整个能谱内积分，计算出探测点总的剂量值。积分公式为

本发明提供了一种用八叉树方法将3dsMax核设施模型体素化后利用点核积分方法计算模型辐射场剂量的仿真方法。

本发明开发出了一种针对于有确定参数的3dsMax核设施模型，使用八叉树方法对模型体素化，选择源项与屏蔽项体素，用点核积分方法计算辐射场剂量的仿真方法。

本发明的有益效果在于：

1、本发明实现了在3dsMax软件内直接对核设施尺寸、材质与能量参数的快速赋值；

2、本发明实现了对复杂模型进行的辐射场剂量计算。

附图说明

图1 3dsMax模型辐射场剂量计算流程图；

图2 3dsMax模型体素化流程图；

图3基于SVR的多层屏蔽累积因子算法模型；

图4点核积分方法的几何模型；

图5辐射场伽马剂量计算流程图。

具体实施方式

下面结合附图举例对本发明做进一步描述：

本发明采用3dsMax软件建模与C++语言编程，主要功能为：构建确定参数的核设施模型，使用八叉树法对模型体素化，最终利用点核积分方法计算辐射场剂量，实现辐射仿真。整个软件包括3dsMax建模、八叉树程序、输入卡的导入、点核积分程序等几个模块。

本发明采用以下技术方案：

1、用3dsMax软件根据确定的核设施参数构建模型，并将文件保存为3DS格式。

在3dsMax软件材质编辑器内直接对材质球命名，然后将材质球导入到使用该材质的实体上。材质的命名规则为

源项类材质命名区分符为SO，命名格式如下：

SO,能量,活度,±份额,其他信息(材料名称等)

举例：SO,0.052,3.27E+06,-1

屏蔽类材质命名区分符为SH，命名规则如下：

SH,密度,原子序数,±份额,其他信息(材料名称等)

举例：SH,2.702,13,-1.0,Al

其中“+”代表原子个数的份额，“-”代表质量份额。

2、导入3DS核设施模型文件，获得模型参数。

定义3DS文件的读入类CLoad3DS，CLoad3DS类包括多个读入函数，如块的读入，顶点数据的读入，面数据的读入，材质的读入等。

读取块。CLoad3DS类的成员函数Import3DS、ReadChunk、ProcessNextChunk用来读取文件并判断其是否是3DS格式，ProcessNextObjectChunk、ReadVertices、ReadVertexIndices、ReadColorChunk、ProcessNextMaterialChunk等成员函数用来读取模型实体、顶点、面、颜色以及材质等内容。函数Import3DS是文件的入口，具体的算法如下：

1)利用fp(fp＝fopen(filename，“rb”))函数以只读方式打开一个3DS文件；

2)验证文件类型，判断是否是3DS文件，若ID是0x4d4d，则为3DS文件，然后调用ReadChunk、ProcessNextChunk等子块读入相应的3DS文件内容，否则返回false；

3)利用fclose(fclose(fp))函数关闭文件。

读取顶点。函数ReadVertices用来读取3DS文件顶点列表的内容，具体的算法如下：

1)定义3个动态指针，分别指向3个动态数组，用来存储模型的所有顶点，定义一个变量count，用来保存模型的顶点总数；

2)将所有的顶点读入到指定的内存；

3)删除指针变量，释放内存。

读取面。函数ReadVertexIndices用来读取3DS文件的面列表内容，具体的算法如下：

1)定义一个动态指针，指向一个动态数组，用来存储模型的面数据，定义一个变量count，用来存储三角形面的数量；

2)将所有的面内容读入到指定的内存；

3)删除指针变量，释放内存。

3、用八叉树法将3dsMax核设施模型分解为体素(图2)。

1)根据模型参数构建包围盒。包围盒为正方体，其边长为整数。设定最小体素边长a＝0.1m，八叉树的最大递归深度为n。初始时将模型所有顶点坐标相加除以顶点个数，获得包围盒中心点坐标。遍历所有顶点，获取与中心点在x、y、z轴方向上的最大距离，并选取其中最大值的2倍为包围盒的拟定边长L。根据条件0.1×2^n-1<L≤0.1×2ⁿ，解得最大递归深度n，并获得包围盒的边长2ⁿ。将包围盒设定为根立方体，递归深度为0。

2)定义体素类COctree，用于记录模型的体素信息。

3)以构成3dsMax模型表面的每个三角形为基本单位，获得每个三角形的顶点和中点坐标。

4)将立方体均分为八等份，并将立方体所装的三角形点元素全部分担给八个子立方体。第一个子立方体的序号i为0。

5)递归深度加1，读入子立方体i。

7)如果i小于8，判断分配到子立方体i的点元素数量是否为零。如果为零，则递归深度减1，子立方体i的标示符为false，i加1，重复步骤5。如果不为零，则判断是否达到最大递归深度。如果达到最大递归深度，则递归深度减1，子立方体标示符为true，i加1，重复步骤5，否则重复步骤4。

8)如果i不小于8，判断递归深度是否大于1。如果递归深度大于1，则递归深度减1，读入下一个子立方体，重复步骤4。如果递归深度小于等于1，程序结束。

4、将确定的体素参数与材质信息写成输入卡。输入卡由6部分组成(表1)，包括问题描述卡、基本体数据卡、材料数据卡、区域卡、放射源卡以及探测点卡，各个卡片之间用回车分割，整个卡片结束用“END”标识。

表1输入文件的数据结构

5、将输入卡导入到点核积分程序内。

6、计算累积因子。

1)单层屏蔽

1991年美国国家标准文件ANSI/ANS-6.4.3给出了23种元素、水、空气与混凝土在光子能量为0.015～15Mev，且光子穿透深度在40mfp以内的累积因子以及相应G-P拟合公式参数。不同材料与不同伽马光子能量的累积因子已经被编制为表格，数据表中不存在的元素可以通过插值计算得到相应的累积因子。对于混合物或者化合物组成的单层屏蔽材料，首先求出等效原子序数，再根据等效原子序数求出G-P公式拟合参数，最后由拟合参数计算出累积因子。在求解等效原子序数过程中，先求出混合物或者化合物的线康普顿衰减系数(σ)与总的线衰减系数(μ)的比值R(σ/μ)，再根据以下公式插值计算出等效原子序数。

式中，Z_eq为等效原子序数；Z₁与Z₂为混合物(化合物)中元素的原子序数；R₁、R₂与R分别为Z₁、Z₂与混合物(化合物)的线康普顿衰减系数与总的线衰减系数的比值。

2)双层屏蔽

计算双层屏蔽结构的累积因子采用经验公式

B＝max{B_a(E,μ_a(d_a+d_b)),B_b(E,μ_b(d_a+d_b))}

式中，能量为E的伽马光子垂直入射到a、b两种不同的屏蔽材料上，若两种材料原子序数相差不大，则可以取两种材料中最大的累积因子作为双层屏蔽总的累积因子，B是总的累积因子；B_a与B_b分别是屏蔽材料a与b的累积因子；μ_a与μ_b分别是屏蔽材料a与b的线衰减系数；d_a与d_b分别是屏蔽材料a与b的厚度；如果两种屏蔽材料的原子序数相差很大，则总的累积因子与屏蔽材料的顺序有关。当低原子序数(Z)介质在前，高Z介质在后，则总的累积因子用高Z介质的累积因子代替，公式为

B＝B_h(E,μ_hd_h)

式中，B_h、μ_h与d_h分别为高Z介质的累积因子、线衰减因子与厚度。当高Z介质在前而低Z介质在后，则总的累积因子与光子能量有关。此时，若光子能量较低，则总的累积因子为两种屏蔽材料各自累积因子的乘积，公式为

B＝B_l(E,μ_ld_l)·B_h(E,μ_hd_h)

若光子能量较高，且超过高Z介质线衰减系数最小值对应的能量(E_h,min)时，累积因子为

B＝B_l(E_h,min,μ_ld_l)·B_h(E,μ_hd_h)

式中，B_l、u_l与d_l分别为低Z介质的累积因子、线衰减因子与厚度；B_h、μ_h与d_h分别为高Z介质的累积因子、线衰减因子与厚度。通常高Z介质的参数E_h,min在3Mev～4Mev之间。

3)多层屏蔽

在计算4层以内(包含4层)屏蔽结构的累积因子，采用布罗杰尔及其合作者根据实验研究的累积因子计算方法，公式如下。

式中，N为屏蔽层数目，B是N层非均匀屏蔽介质的累积因子，B_n是第n层屏蔽材料组成的均匀介质的累积因子。

对于大于4层屏蔽结构的累积因子采用Suteau与Trontl提出的迭代算法计算。首先，引用数据库，包括单层、双层以及多层屏蔽情况。然后，确定输入空间，包括屏蔽层前n层(n＝1,2,…,N)的有效累积因子与第N层屏蔽材料的累积因子、屏蔽层前n层(n＝1,2,…,N-1)的有效原子序数与第N层屏蔽材料累积因子、屏蔽层前n层(n＝1,2,…,N-1)的等效屏蔽厚度与第N层的等效屏蔽厚度、光子能量。其中，等效屏蔽厚度为前n层屏蔽厚度之和，有效原子序数采用如下公式计算。

式中，Z_eff是等效屏蔽层原子序数，Z₁/Z₂与d₁/d₂分别是屏蔽材料原子序数与厚度。最后，Trontl以支持向量机器算法库LIBSVM为计算工具，并采用迭代模型(图3)实现多层屏蔽累积因子计算。

7、计算伽马射线在辐射场中的平均自由程。

平均自由程计算公式为

式中，t(E)是伽马光子从点源到探测点穿过所有屏蔽材料的平均自由程，i为伽马射线穿过的空间区域编号；ρ_i为空间区域i的材料密度；μ_i(E)/ρ为在光子能量为E时，空间区域i的材料的质量减弱系数；d_i为伽马射线在区域i中的几何距离。

计算平均自由程的关键是求出射线穿过的屏蔽材料以及在各屏蔽层(区域)中穿过的几何距离。本发明采取以下步骤实现射线穿过的屏蔽层以及平均自由程的计算。

1)计算点核到探测点的伽马射线与所有基本几何体的交点，并去掉交点中的重复项。在计算过程中，若点核与探测点都在某基本几何体内部，则伽马射线与该基本几何体没有交点，此时直接计算射线与下一个基本几何体的交点，利用这一过程可以显著地提高计算效率。

2)根据交点与点核的距离，将所有交点按照由近到远的顺序进行排列。若交点数目少于2则不用排序。

3)将所有交点进行排序后，在交点数据结构的首末位置分别添加点核与探测点的坐标，然后计算相邻两个交点的中点坐标以及相邻两个交点的距离。

4)依次判断中点所在的区域，根据区域的材料确定屏蔽层，与中点对应的距离即是屏蔽层的厚度。若相邻两个屏蔽层的材料相同，则合并该屏蔽层。

5)根据伽马光子能量求出各屏蔽材料的线减弱系数，并利用平均自由程计算公式求出伽马射线平均自由程。

8、利用SQLite数据库引擎建立通量率-剂量率转换因子、各化学元素与材料的质量衰减因子以及单层累积因子数据库。

通量率-剂量率转换因子采用ANSI/ANS-1977与ICRP-21报告中提供的数据创建Convertor.db数据库。在光子能量在0.015～15Mev区间内，利用WinXCom程序计算化学周期表前100号元素、空气、水以及标准混凝土的质量衰减因子，并以此创建Mafactor.db数据库。单层屏蔽累积因子以美国国家标准文件ANSI/ANS-6.4.3-1991提供的数据为基础，得到23种元素、水、空气与混凝土在光子能量为0.015～15Mev内的累积因子数据库Bufactor.db。

9、运用布尔连接运算符对箱体进行组合运算，构造复杂的辐射场几何结构。

在伽马剂量计算模块中，利用布尔运算符(并、交与补)将箱体组成复杂的区域来描述辐射场。本发明中必须给出箱体(体素)的某一个顶点坐标P(x₀,y₀,z₀)，以及能够表示箱体长宽高且相互垂直的三个向量V1(v_1x,v_1y,v_1z)、V2(v_2x,v_2y,v_2z)与V3(v_3x,v_3y,v_3z)，才能够唯一确定箱体的位置与尺寸。组合几何方法的主要过程如下所述。

1)创建基本体。根据输入卡中箱体(体素)在空间中的位置与尺寸，按照编号的先后顺序建立各箱体。

2)创建子区域。根据输入卡内箱体的“+”或者“-”，构成一个子区域。其中“-”表示在箱体内部，“+”表示在箱体外部。

3)创建区域。若干子区域用布尔逻辑符“AND”与“OR”连接构成区域。其中“AND”表示区域是两个子区域的交集(通常省略AND)，“OR”表示区域是子区域的并集。至此，可以描述复杂的辐射场几何结构。

10、用点核积分方法计算三维辐射场剂量(图5)。

1)导入利用组合几何方法构造的输入卡。

2)计算放射源各能量点的通量率-剂量率转换因子，各屏蔽材料的线衰减因子及G-P拟合参数。

3)将放射源离散为若干点核。

4)计算点核到探测点的平均自由程及屏蔽层。

5)若点核与探测点之间只有一层屏蔽，则采用G-P拟合公式计算累积因子，若屏蔽层数目在4层以内(包含4层)采用经验公式计算累积因子，若屏蔽层数目大于4则采用迭代模型。

6)计算累积因子并根据点核公式计算点核在探测点的剂量率。

7)完成所有点核的计算，并将所有点核在探测点的剂量率求和。

8)完成所有放射源计算，并将所有放射源在探测点的剂量率求和。三维辐射场剂量计算完成。

Claims

1.一种基于3dsMax的核设施模型辐射场剂量仿真方法，其特征是：

(2)导入3DS核设施模型文件，获得模型参数；

(3)用八叉树法将3dsMax核设施模型分解为体素，具体包括：

1)根据模型参数构建包围盒，包围盒为正方体，其边长为整数，设定最小体素边长a＝0.1m，八叉树的最大递归深度为n，初始时将模型所有顶点坐标相加除以顶点个数，获得包围盒中心点坐标，遍历所有顶点，获取与中心点在x、y、z轴方向上的最大距离，并选取其中最大值的2倍为包围盒的拟定边长L，根据条件0.1×2^n-1<L≤0.1×2ⁿ，解得最大递归深度n，并获得包围盒的边长2ⁿ，将包围盒设定为根立方体，递归深度为0；

2)定义体素类COctree，用于记录模型的体素信息；

3)以构成3dsMax模型表面的每个三角形为基本单位，获得每个三角形的顶点和中点坐标；

4)将立方体均分为八等份，并将立方体所装的三角形点元素全部分担给八个子立方体，第一个子立方体的序号i为0；

5)递归深度加1，读入子立方体i；

7)如果i小于8，判断分配到子立方体i的点元素数量是否为零，如果为零，则递归深度减1，子立方体i的标示符为false，i加1，重复步骤5)；如果不为零，则判断是否达到最大递归深度，如果达到最大递归深度，则递归深度减1，子立方体标示符为true，i加1，重复步骤5)，否则重复步骤4)；

8)如果i不小于8，判断递归深度是否大于1，如果递归深度大于1，则递归深度减1，读入下一个子立方体，重复步骤4)；如果递归深度小于等于1，结束；

(4)将确定的体素参数与材质信息写成输入卡；

输入卡由6部分组成，包括曲面卡、基本体数据卡、材料数据卡、区域卡、放射源卡以及探测点卡；

(5)将输入卡导入到点核积分程序内；

(6)计算累积因子；

(7)计算伽马射线在辐射场中的平均自由程；

(10)用点核积分方法计算三维辐射场剂量，具体包括：

设核反应堆以及其他放射源由一系列离散的各向同性点源组成，并且通过对各离散点核的贡献求和来获取放射源在探测点总的辐射效应，在点核积分方法的几何模型中，各点核在探测点的剂量值为

式中，r_p与r_d分别为点核与探测点的位置；E为光子能量；C(E)为伽马光子辐射效应转换因子；S(E)是点核源项强度；B(E,t(E))是累积因子；t(E)是伽马光子从点源到探测点穿过所有屏蔽材料的平均自由程，计算公式为

式中，i为伽马射线穿过的空间区域编号；ρ_i为空间区域i的材料密度；μ_i(E)/ρ为在光子能量为E时，空间区域i的材料的质量减弱系数；d_i为伽马射线在区域i中的几何距离；

将探测点的剂量值在整个源项体积与整个能谱内积分，计算出探测点总的剂量值；积分公式为

2.根据权利要求1所述的基于3dsMax的核设施模型辐射场剂量仿真方法，其特征是：3DS核设施模型文件由许多块组成，每个块首先描述其信息类别；块的前两项信息分别是：块的ID和块的长度，块的ID作为块的标示，块的长度表示下一个块相对于该块起始位置的偏移字节数；3DS文件最开始出现的主块是基本块，它包含了整个文件；3D编辑程序块ED—IT3DS定义物体的形体数据，块中包含了一个物体描述子块EDIT_OBJECT，其下有三角形列表子块OBJ_TRIMESH；OBJ_TRIMESH包括了体素化程序所需的模型表面信息。

3.根据权利要求1或2所述的基于3dsMax的核设施模型辐射场剂量仿真方法，其特征是所述计算累积因子具体包括：采用ANSI/ANS-6.4.3数据库以及G-P拟合公式计算单层累积因子，采用经验公式计算双层或多层材料累积因子。