中子输运计算集成软件平台的研究与设计

李 磊, 王世庆, 李 伟, 柳 建

(1.核工业西南物理研究院, 四川 成都 610225;2.中广核研究院有限公司, 广东 深圳 518000)

在中子输运分析中，主要的方法有确定论方法以及非确定论方法[1]，其中非确定论方法即随机模拟方法,具有适用复杂几何机构、能进行精细模拟等优势，随着现代计算机计算能力的大幅度提升，在核工程领域得到越来越多的应用[2]。由美国阿拉莫斯实验室(LANL)开发的蒙特卡罗程序MCNP是随机模拟程序中的佼佼者，已广泛应用于辐射防护、反应堆设计、临界装置实验等领域[3]。但是MCNP输入卡使用文本结构，其内容全部需要人工输入，效率低并且容易出现输入卡编写错误，对工程人员的经验要求较高，在复杂几何模型中，这些问题更加凸显。

为解决MCNP输入卡复杂、效率低、易出错的问题，一些学者对MCNP的建模自动化进行了研究，并且开发了特定功能的应用工具。从最早的Visual Editor、Sabrina和Moritz[4]等模型可视化工具到MCAM[5]，极大地提高了程序使用的便利性。这些工具使用的方法也有多种，例如罗月童等人提出的BREP到半空间转化算法[6]，张建生等人提出的UG模型空间树转化算法[7]，王寒冰提出的基于特征的BREP到CSG模型转换算法[8]以及吴烔等人运用卷积神经网络(CNN)对CAonSTEP算法进行改进的ICAonSTEPS算法[9]等。

上面提到MCNP辅助工具都是基于大型三维建模软件的二次开发，属于重量级工具。随着“第四代核能系统”研究不断深入，需要对不同堆芯模型进行设计验证，为满足核反堆数字化设计需求[10]，将更多工作交由计算机完成。因此，亟需开发一套集成MCNP输入卡自动生成、自动调用计算核心程序及计算结果后处理等功能于一体的轻量级软件平台，以降低MCNP程序使用门槛和时间成本，提高工作效率。本文工作依托于国家重点研发计划“核安全与先进核能技术”重点专项，基于堆芯模型快速搭建、MCNP灵活输入以及计算结果提取分析等功能，实现反应堆堆芯物理特性的快速分析[11]，提高项目研发效率。

1 集成软件平台的功能及关键模块

集成软件平台主要由三部分组成，包括输入卡自动生成模块、计算核心调用模块以及结果处理模块。在整个中子输运计算过程中，前后处理占了60%～80%的时间[12]，包括几何模型构建、计算参数设置以及结果可视化分析等，相比于传统的单独处理方法，将其功能模块集成后能节省大量时间。输入卡生成模块将实体几何模型转化为MCNP输入卡中的几何模型信息，再补全剩余卡片信息生成完整输入卡文件。计算核心调用模块中将上一模块生成或外部空间导入的输入卡文件使用MCNP程序进行计算，得到计算结果文件。结果处理模块中进行计算结果的数据处理和分析，其处理的结果文件可以是计算模块中产生的，也可以从外部空间中导入。整个集成软件平台的功能流程图如图1所示。下面对其中的几个关键技术进行阐述。

图1 集成软件平台的功能流程框图Fig.1 The function flow chart of the integrated software platform

1.1 BREP文件到MCNP几何输入卡的自动生成

经过对多种实体几何模型文件类型的比较，最终选择了可读性高的边界表示法(BREP)文件作为实体几何模型源文件。BREP文件使用边界来表示实体几何模型，通过基本几何元素(点、线、面、体等)来存储几何信息，同时依据已知的拓扑关系(体→面→环→边→点)来构建各基本几何元素间的连接关系，进而实现对实体模型的表示。BREP格式文件以文本形式存储，便于读写。

自动生成输入卡模块程序将BREP文件转化为MCNP输入几何卡，然后添加材料卡，数据卡等信息，生成完整的MCNP输入卡文件。对如图2所示的几何模型，其BREP文件部分内容如图3所示。

从图3中可以看到，BREP文件中包含了实体模型的几何信息和拓扑信息。几何信息中位置信息部分以标识“Locations”开始，其后的数字表示位置信息的数量。每个位置信息是一个 3×4的矩阵，描述三维空间的线性变换。“Surfaces”标识符下包含几何模型所有的曲面信息，图3中的曲面信息为一个圆柱面的表示方式，三维正交坐标系中圆柱面的轴通过点(-6,-6,0)，方向为[0,0,1]，半径为0.5，其参数方程为：

S(u,v)=P+r·(cos(u)·Dx+sin(u)·Dy)
+v·Dv,(u,v)∈[0,2π)×(-∞,∞)

(1)

其中，Dv，Dx，Dy一起组成三维正交坐标系，圆柱面的中心轴通过点P，方向为Dv，圆柱面的半径为r。“TShapes”标识符下包含几何模型中各基本几何元素间的拓扑关系信息，因此可以获取几何实体与各曲面之间的关系。“Fa”表示面(face)，“So”表示实体(solid)，图3中的拓扑信息表示一个由圆柱面和上下底面构成的圆柱体。

MCNP输入卡中的几何描述包括曲面卡和栅元卡两部分。曲面卡包含几何模型的所有曲面信息，相当于BREP文件中的“Surfaces”部分，因此直接将曲面信息转化为MCNP输入卡中的栅元卡信息，MCNP中曲面的表示规则见表1， MCNP中使用一般方程对曲面进行描述，而BREP中对曲面的描述则使用参数方程，因此生成MCNP的曲面卡信息需要获得两者之间的转换关系。由于球面、圆柱面等曲面与坐标轴之间存在关系，其一般方程参数值可以根据BREP文件中的曲面参数可以直接确定，本文主要处理平面的参数方程与一般方程之间的关系，根据向量共面的条件可知：

(2)

其中，(x0,y0,z0)为平面上的点，(X1,Y1,Z1)和(X2,Y2,Z2)分别为平面上的向量Du和Dv的坐标。求解行列式再与平面一般方程比较，获得两者之间的转换关系(式3)。

(3)

表1 MCNP曲面描述

栅元卡中的栅元由曲面定义的半空间通过正则运算组合而成，模型中的每个区域都必须定义，不能存在空隙。因此MCNP栅元构建的算法步骤为：

1)根据BREP文件中提取的体与面拓扑信息生成组成几何模型的实体集，若某个实体存在位置变换信息，则将当前位置与位置变换矩阵Q相乘以确定实体的最终位置；

2)确定实体集中各个实体之间的关系，本文中只考虑一个实体完全包含另一个实体和两个实体分开这两种情况，最终生成实体树表示各个实体之间的关系；

3)最后从叶子节点开始遍历实体树上的所有实体节点，生成“a ±f1±f2… #b1#b2…”格式的栅元信息，其中a为栅元号，f1、f2为当前组成实体曲面边界的曲面号，±表示曲面方向，b1、b2为实体节点的所有子实体节点的栅元号。

使用模型转换算法将图2所示的实体模型的BREP格式表示转换为MCNP输入卡的几何描述格式，最终结果如图4所示，符合MCNP输入卡的格式要求。

图4 BREP→MCNP模型转换结果Fig.4 The result of converting BREP file into MCNP geometry input card

1.2 对MCNP程序计算核心的调用

对自动生成的几何卡补充材料、数据等信息后生成完整的输入文件，在集成平台内部调用MCNP计算核心就能够直接进行中子输运计算，并捕获程序输出的计算过程信息，同步显示到当前平台的信息区，运行界面如图5所示。

图5 MCNP计算信息图Fig.5 The MCNP calculation information

1.3 计算结果后处理

MCNP的计算结果是如图6所示的文本格式，包含计算结果数据以及一些特殊字符串[13]。当前对MCNP计算结果的处理绝大多数仍使用传统的手工数据分析方法，需要人员从结果文件中提取数据，再导入专业的数据分析软件中进行处理，工作量大且效率低[14]。集成平台中的后处理模块较好地实现了数据提取和绘制图表两个功能。

图6 MCNP结果文件片段Fig.6 The result file fragment of MCNP

数据提取功能通过检索结果文件中关键字定位数据所在位置来实现，如图6中的“k(coll)”字符串，其后为具体的数据值。但从图6中可以看到，数据值之间可能存在“|”等特殊字符，因此每获取一行数据，需要使用正则表达式过滤数据中的特殊字符并将行数据分隔为一维数组，最终处理完所有数据后获得一个二维数组，将二维数组保存到Excel格式的中间文件中[15]，完成数据提取工作。

绘制图表功能利用Java绘图工具包JFreeChart[16]将结果文件中提取的数据绘制为折线图。折线图绘制过程中，先获取数据集对象DataSet，将绘图数据添加到数据集中，然后调用JFreeChart的API生成并显示折线图表。

2 集成软件平台测试结果

完整的集成软件平台程序的界面如图7所示，左侧为程序中各模块的操作区，各个模块既可以一起使用，也可以单独使用。右侧为信息区，主要显示各模块中生成的各种文件信息、过程信息等。

图7 程序界面Fig.7 The program interface

为测试软件中各模块的功能以及效率，对某小型堆堆芯的临界系数Keff进行计算，以获取最佳燃料富集度。所选取堆芯得具体参数如表2所示[17]，该堆型采用7×7的堆芯布局，包括37个燃料组件，每个燃料组件为标准的压水堆燃料组件[18]，如图8所示，其中B为可燃毒物棒，G为导向管，I为仪表管，余下的为UO2燃料棒。计算中使用富集度为2.4%的UO2燃料棒。

表2 小型堆堆芯几何参数

图8 燃料组件布局Fig.8 The arrangement of the fuel assembly

本计算中子源位于堆芯中心，源强为20 000，模拟300代中子循环下堆芯临界系数Keff的变化，计算完成后提取结果文件中各中子代循环下的Keff值，并保存到Excel文件中，提取数据如图9所示。然后使用绘图功能生成图10所示的结果曲线图。当程序计算到100代循环时，临界系数Keff已趋于稳定，因此为节省时间，后续计算取100代中子循环的结果为有效结果进行分析。

图9 临界计算数据Fig.9 Results of criticality calculation

图10 堆芯临界计算结果Fig.10 Calculation results of core criticality

可以看出，本软件处理数据较为方便，并且后续容易根据用户需求开发更丰富的后处理功能。

从图10可以看出，当前计算模型下，堆芯临界系数在1.07上下波动，此时堆芯处于超临界状态。为了获得临界状态，在相同几何布局下，调整燃料富集度分别为1.0%、1.5%、2.0%，进行了多轮计算，以获得堆芯达到临界状态时的燃料富集度。对比计算结果如图11所示。

图11 不同富集度下的临界系数Keff比较Fig.11 Comparison of the critical coefficient Keff under different enrichment degrees

由图11可知，随着燃料富集度的提升，堆芯临界系数Keff逐渐增大，该堆型的堆芯达到临界状态时的燃料富集度在1.5%～2.0%。在此富集度区间再进行多轮计算，最终得到该堆型堆芯临界的最佳燃料富集度为1.85%，结果如图12所示。

图12 富集度1.85%的Keff变化曲线Fig.12 Variation of the Keff in every generation with the enrichment of 1.85%

整个测试过程显示，采用集成软件平台，完成完整一轮堆型验证比原有方式节约2人·天，总体效率提高约60%。

3 结束语

通过对BREP文件格式和MCNP程序输入输出文件的深入研究，采用Python和Java联合编程技术，将几何模型转输入卡的处理、中子输运计算核心的调用、计算结果的图形化后处理实现了一体化集成，搭建了中子输运仿真集成软件平台，平台中各个模块既可以关联使用，又能单独使用，具有较高的灵活性。

通过堆芯临界分析项目对集成软件平台的有效性进行了验证。结果显示，该软件平台在MCNP输入文件准备中，解决了人工建模易出错、耗费时间长的问题，在后处理中解决了数据处理繁琐的问题，整个仿真流程的效率提高60%，并且降低了对工程分析人员的技术要求。这为满足日益增加的反应堆物理随机分析场景，提供了一个便捷有效的工具。