球谐函数有限元程序NECP-FISH的开发及其在聚变堆包层中子学分析中的应用

苗建新，曹良志,*，方 超，贺清明，曹启祥，尹 苗

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.核工业西南物理研究院，四川 成都 610041)

采用球谐函数展开中子输运方程的角度变量，可避免离散纵标方法中出现的射线效应；采用有限元方法离散输运方程的空间变量，可减少传统结构网格对复杂几何建模上的近似；将两种方法结合起来求解中子输运方程可兼具二者的优点。标准的Galerkin有限元方法适用于椭圆方程的求解，如二阶中子输运方程，扩散方程等；对于一阶中子输运方程这种双曲方程，直接采用Galerkin有限元方法会出现非物理的震荡现象。国内外针对球谐函数有限元方法的研究主要求解的是二阶自共轭方程[1]和二阶偶对称方程[2]，但由于截面会出现在方程分母上，对于包含真空区域的问题无法进行求解，限制了方法的应用场景。Buchan等[3]针对求解一阶中子输运方程的球谐函数有限元方法进行了研究，解决了数值不稳定现象，但目前还未应用于三维大规模问题。

包层作为聚变堆中的重要结构，承担着氚增殖、能量转换和辐射屏蔽等一系列功能。因此对包层的中子学分析显得尤为重要。聚变堆包层由第一壁、上下盖板、增殖单元、冷却剂管道以及背板等结构组成，存在较多的非规则几何结构，如采用基于结构网格的确定论方法进行计算要引入大量的近似，计算精度较低。蒙特卡罗方法虽能描述复杂的几何结构，但计算耗时较长。由于有限元方法对于复杂几何的处理具备天然的优势，近年来，基于有限元方法的确定论程序开始应用于包层中子学分析。如美国的商用软件ATTILA[4]采用离散纵标和间断有限元方法，已应用于国际热核聚变实验堆(ITER)的中子学分析中。韩国的AETIUS[5]与ATTILA采用相同的方法，已应用于韩国氦冷陶瓷试验包层的中子学分析中。目前国内进行聚变堆包层中子学计算大多采用蒙特卡罗程序，还未开发能应用于包层中子学分析的有限元程序。

本文研究求解一阶中子输运方程的球谐函数稳定有限元方法，并基于此方法开发中子学分析程序NECP-FISH。最终将其应用于氦冷陶瓷包层，并与蒙特卡罗程序NECP-MCX[6]的计算结果进行对比分析。

1 球谐函数有限元方法

NECP-FISH程序采用球谐函数和多尺度稳定有限元方法求解一阶稳态中子输运方程。

不考虑裂变材料，中子输运方程右端只有外源项和散射源项，多群输运方程不同能群通过散射源项进行耦合，此时方程可写为单能形式：

(1)

(2)

根据加权残差法获得一阶输运方程的弱形式：

〈w,φ〉∂V+(L*w,φ)V=(w,s)V

(3)

本文所采用的多尺度有限元方法是子网格有限元方法中的一种。子网格有限元方法最早由Hughes[7]提出，Buchan在2008年第1次将其应用到求解中子输运方程中[3]。此方法最重要的思想就是将未知量分解为两部分：

φ=φc+φr

(4)

式中：第1部分φc为连续可解部分，定义在全局尺度上；第2部分φr为残差部分，这部分在最终离散的方程中并不进行求解。根据上述对未知量的分解，同时由于残差部分定义在每个单元上，边界条件只应用于连续部分，所以式(3)的弱形式可写为：

〈w,φc〉∂V+(L*w,φc+φr)=(w,s)

(5)

同时残差方程可表示为：

Lφr=s-Lφc

(6)

此残差方程可近似求解，使用不同的近似方法产生了不同的子网格有限元方法。本文中的代数子网格方法(algebraic sub-grid scale method)采用代数方程对通量密度的残差部分进行近似：

φr=λ(s-Lφc)

(7)

式中，λ相当于对L-1的近似。将式(7)代入到式(5)中，可获得代数子网格方法的弱形式为：

〈φ,φc〉+(L*φ,φc)-(L*φ,λLφc)=

(φ,s)-(L*φ,λs)

(8)

上式中的未知量φc可由基函数展开φ(r,Ω)=φT(r,Ω)φ，此基函数可用分段多项式函数N(r)与球谐函数Y(Ω)的克罗内克积的形式来表示：

φ(r,Ω)=N(r)⊗Y(Ω)

(9)

最终方程的离散形式可写为：

(10)

在每个内部单元和边界单元上，式(10)中左端的每项均可计算得到并组装到整体刚度矩阵中，最终通过求解线性方程组来获得每个网格节点的通量密度。对此球谐函数有限元方法的详细推导见参考文献[8]。

2 程序开发及功能介绍

基于上述理论方法，并借助开源的前后处理平台SALOME[9]，开发了中子学分析程序NECP-FISH，程序的整体框架如图1所示，该程序具备从可视化建模、中子学分析到结果可视化展示的一系列功能。其中中子学分析涵盖了中子通量密度、氚增殖比、中子释热的计算，主要由球谐函数有限元求解器实现。借助SALOME平台开发程序的优点主要是：首先，SALOME平台集成了很多开源工具(如实现CAD建模的OpenCASCADE，网格剖分的Gmsh和NETGEN，可视化的ParaView)，借助其开发程序，可使程序的功能更完备；其次，SALOME平台提供了用户二次开发的接口，用户通过开发接口程序，可方便地将自己的数值计算程序接入到平台中形成完整的分析程序，可实现数据的内部传递，方便用户操作；最后，SALOME平台的使用也较为简单，它在Windows和Linux系统下均可使用。借助其开发的程序兼容性也较好。但另一方面，由于平台集成的多为开源工具，相比于提供前后处理功能的商用软件，在功能上还有一定欠缺。

图1 NECP-FISH程序框架

2.1 可视化建模

程序的可视化建模通过SALOME平台集成的CAD和网格剖分工具实现。程序可直接导入设计好的几何文件进行建模，也可对分开导入的几何模型进行装配从而建立最终的计算模型。对于结构简单的基准题，用户还可直接构建几何体再通过布尔运算进行建模。对于建立好的几何模型，NECP-FISH程序可采用不同的网格剖分方法，目前程序有限元求解器支持的非结构网格类型为二维三角形和三维四面体网格。利用程序建立的CFETR氦冷陶瓷包层的中子学分析模型如图2所示。

图2 NECP-FISH程序建立的模型

上述的建模和剖分网格操作均可通过程序提供的用户图形界面可视化完成，同时，为更方便完成中子学计算，NECP-FISH程序在SALOME平台中采用Python结合Qt库的方式开发了接口程序，通过接口程序的可视化界面，将多群截面、求解参数、边界条件、中子源以及剖分出的网格存储在HDF5格式的输入文件中。相比于传统程序的文本输入，HDF5文件的特点是可分级分类存储数据，其存储相同大小的数据量所占的内存空间更小。最终生成的输入卡片可直接提供给有限元求解器进行计算。

2.2 中子学分析

程序的中子学分析功能是由球谐函数有限元求解器(PN-FEM)实现的，有限元求解器可计算出各网格节点上的中子通量密度，再根据氚增殖比以及释热率的计算公式可给出相应物理量的计算结果。

1) 有限元求解器

NECP-FISH程序的核心是有限元求解器。求解器是基于函数库进行开发的，图3所示为求解器的整体框架。Fortran通用函数库和输入输出模块将输入文件的内容读入，再通过球谐函数模块、有限元模块、材料截面和源模块构造出要求解的问题。问题的求解主要是由开源数学库psblas[10]组成的线性代数求解模块实现的。通过建立不同模块的函数库，可使程序结构清晰易于理解，同时基于函数库开发程序也更快捷。

图3 求解器的整体框架

求解器的中子学计算流程如图4所示，整体可分为预处理、求解和后处理。在预处理过程中读取输入并根据计算核数决定是否进行空间区域分解。求解过程主要是进行多群迭代，多群迭代是能量由高到低依次对每个能群进行求解。由于计算的问题是屏蔽问题，所以不存在裂变源迭代。同时由于群内散射移到了方程左端，所以求解过程中不需要散射源迭代。程序利用前面能群计算出的通量密度矩更新散射源后即可进入单群方程求解，之后在每个单元上计算单元矩阵和向量再组装到整体刚度矩阵和向量中。组装好的线性代数方程组采用广义极小残余算法(GMRES)进行迭代求解。对于没有上散射的问题，只需1次能群扫描，对于包含上散射的问题，需从存在上散射的最高能群进行多次能群扫描，直至多群通量密度收敛。最后在后处理模块中实现氚增殖比、中光子释热等物理量的计算，并生成可视化文件。

图4 求解器的中子学计算流程图

2) 程序空间并行

随着网格数的增多，计算规模迅速增大，因此需利用高性能计算平台并行计算以提高计算效率。由于球谐函数很难实现各角度矩的解耦，所以程序只采用空间并行。并行策略是基于消息传递接口(MPI)的分布式并行。并行求解示意图如图5所示。首先在预处理过程中进行空间区域分解，非结构网格的空间区域分解算法十分复杂，通用的做法是使用专门的程序包。NECP-FISH程序采用的是METIS程序包，将非结构网格各节点按特定编号排列好，再调用相应函数即完成区域分解。这样就将整个计算问题划分为几个小的区域分配到各核上，并明确了核与核之间要通信的网格节点编号。区域分解后相当于在每个计算核心上组装并求解整体线性方程组的一部分，迭代求解过程中的矩阵和向量相乘以及向量点积计算要进行核间的通信。最终统一区域交界处各点的计算结果即完成并行计算。

图5 程序并行求解示意图

2.3 结果可视化

结果的可视化需生成存储计算结果的特定格式文件，该文件包括网格信息、离散节点或单元上的计算结果。NECP-FISH程序可视化结果包括中子通量密度分布以及释热分布，采用XDMF格式。该格式采用XML格式存储数据路径，再读取HDF5文件中的数据来完成结果的可视化。相比于传统的VTK格式文件，XDMF格式的文件可并行输出，同时数据存储所占的内存更小。文件可导入到SALOME平台中的ParaVis模块中查看，也可采用其他可视化软件进行查看。

3 氦冷陶瓷包层计算结果

将NECP-FISH程序应用于聚变示范堆氦冷陶瓷包层的中子学分析中，此包层由核工业西南物理研究院设计[11]。由于包层内部为增殖单元对称排布，且增殖单元涵盖了整个包层的全部材料，因此分析单个增殖单元即可反映整个包层的中子学特性。本文选取具有代表性的包层增殖单元进行计算。图6为增殖单元的材料布置和边界条件，氚增殖剂为Li4SiO4，中子倍增剂为Be，两者做成U型球床的结构，中间用冷却板间隔开。背板及冷却板的材料均采用低活化的铁素体钢(CLF-1)。在第一壁前增加了2 mm的钨涂层，同时在钨涂层前定义了各向同性的均匀体源。

根据图6的几何布置和尺寸，利用NECP-FISH程序建立计算的三维模型，并剖分出非结构的四面体单元，如图7所示。

图6 增殖单元的材料布置和边界条件

a——建立的三维增殖单元模型；b——程序剖分的网格

分别采用NECP-FISH和NECP-MCX对此模型进行计算。NECP-FISH的多群数据库采用FENDL-2.1的MATXS格式数据库[12]，并通过TRANSX处理获得多群宏观截面，能群结构为Vitamin-J 175群的能群结构。散射截面展开阶数为P4，球谐函数展开阶数为P15，剖分的四面体网格数为292 593。NECP-MCX采用FENDL-2.1的ACE格式数据库，投入了100组，每组1 000 000个粒子。图8为NECP-FISH计算出的中子通量密度的空间分布，随着与中子源距离的增大，中子通量密度下降了2个数量级。

图8 NECP-FISH计算获得的中子通量密度空间分布

各区域的平均积分通量密度与NECP-MCX计算结果的对比如图9所示，通量密度的偏差在距源较近的区域稍大，所有区域的相对偏差均在6%以内。同时，图中还给出了蒙特卡罗程序MCNP[13]的计算结果，两蒙特卡罗程序的相对偏差均在0.4%以下。

图9 各区域平均积分通量密度

对此问题进行氚增殖比的计算，氚增殖比的定义为1个归一的聚变中子在包层中产生的氚原子个数，产氚主要是通过中子与包层中的含锂材料发生反应实现的，因此NECP-FISH程序计算氚增殖比的公式为：

(11)

式中：Σt,Li4SiO4(E)为硅酸锂材料的宏观产氚截面；Sn为中子源强。NECP-FISH程序计算的此增殖单元的整体氚增殖比为1.339，NECP-MCX的计算结果为1.332，相对统计方差为0.006 7%。MCNP的计算结果为1.334。氚增殖比的分布及NECP-FISH与NECP-MCX程序计算结果的相对偏差列于表1。

表1 氚增殖比的计算结果

各区域中子释热的结果列于表2，中子释热是由通量密度与截面相乘再进行积分获得的：

表2 增殖单元各区域中子释热对比结果

Kn=∬Σh(E)φ(r,E)drdE

(12)

式中，Σh(E)为宏观释热截面。由上式可知，中子释热偏差与通量密度的偏差直接相关，NECP-FISH与NECP-MCX计算结果的最大偏差出现在靠近中子源的区域，所有区域的偏差都在6%以内。在包层核热计算方面，前文所提到的确定论程序ATTILA与蒙特卡罗程序最大相对偏差达到了28%[14]。

4 结论

本文研究了球谐函数有限元方法，并基于此方法开发了中子学分析程序NECP-FISH。该程序具备可视化建模、中子学计算和结果可视化功能。利用NECP-FISH对聚变堆氦冷陶瓷包层进行了建模，并计算了中子通量密度、氚增殖比以及中子释热3个重要的中子学参数。结果表明，NECP-FISH与蒙特卡罗程序NECP-MCX计算的TBR的相对偏差为0.56%，各区域平均积分通量密度最大相对偏差为钨涂层的-5.27%，其余区域均在5%以内，中子释热最大相对偏差为-5.79%。结果证明了NECP-FISH程序能应用于聚变堆包层中子学分析。某些区域偏差较大的原因可能是使用的多群数据库与连续能量数据库的差别，后续将进行进一步的研究。