基于二十面体间断有限元离散求积组的IRI-TUB基准题验证

代 妮，张 斌，陈义学

(华北电力大学 核科学与工程学院，北京 102206)

粒子输运计算在核装置的屏蔽设计中起到至关重要的作用，屏蔽计算的核心是粒子输运方程的求解。稳态输运方程是一个包含能量、空间和角度共6个自变量的线性微分-积分方程，通常采用数值离散的方法近似求解[1-2]。离散纵标法(SN)作为应用最广泛的确定论输运计算方法之一，具有计算速度较快、精度较高、适合求解深穿透输运问题等优点[3]。SN对输运方程的角度变量直接离散，将原来连续的方向转化为有限个离散方向进行求解。对于各向异性较强的深穿透屏蔽问题，通量密度的角度分布不光滑甚至出现间断，传统求积组无法准确对高阶球谐函数积分，导致离散误差较大，严重影响了屏蔽计算的精度与可靠性。因此，研究具备积分分段间断函数能力的求积组是屏蔽计算中的一个重要研究方向。

在SN的求积组研究领域，1964年Carlson和Lathrop[4]提出层对称求积组(LS)。该求积组按层划分，层的交点即为求积点，根据矩条件求解方向余弦与权重系数。LS是目前应用最广泛的求积组，但当离散阶数大于20时，其权重系数会为负。随后Carlson[5]又提出等权重求积组(EQN)，人为地使各离散方向的权重系数相等以避免权重为负，但不便于构造高阶求积组。1987年Walters[6]提出勒让德-切比雪夫求积组(PNTN)。该求积组具有积分精度高、易于构造大量离散方向等优势。1995年Thurgood等[7]提出全对称TN求积组，采用八面体向球面投影生成求积点，对应球面三角形面积为权重系数，相比于LS精度更高且保证权重恒为正。2007年，Endo和Yamamoto[8]研究了偶阶奇阶求积组(EON)，使用单动量偶阶矩和奇阶矩求解求积组参数。该求积组能同时满足更多矩条件，但仍存在阶数限制，大于16阶会出现负权重。2010年Jarrell等[9]开发了线性间断有限元求积组，根据八面体三角形几何特性定义离散方向，通过对线性间断有限元基函数积分得到相应权重。该求积组可在1/8球面上准确积分0阶和1阶球谐函数，能产生大量离散方向且易于实现局部细化。2012年Ahrens[10]根据旋转不变性提出了二十面体求积组，该求积组能准确积分高阶球谐函数，但由于其不具有八面体对称性而难以处理反射边界条件。2015年Lau和Adams[11]基于球面四边形研究了线性和二次间断有限元方法及相适应的局部细化求积组。

本文结合间断有限元思想，构造二十面体线性和二次间断有限元求积组，以适应角通量密度各向异性较强的屏蔽问题，减小角度离散误差，提高屏蔽计算方法的计算精度与可靠性。

1 理论模型

1.1 间断有限元求积组构造流程

在离散纵标法中，离散方向及其对应权重系数的集合称为离散求积组。一般来说，采用离散方向的极角θ及方位角γ的方向余弦(μ,η,ξ)定义坐标系。间断有限元求积组的构造流程如图1所示。首先将正二十面体嵌入单位球内，其每个面分别对应1个角度区域。在每个平面三角形内建立坐标系(u,v)并定义初始离散方向，将平面上的离散点投影到单位球表面，即可得到相应且唯一的求积点，如图2所示。其次根据初始离散方向求解间断有限元基函数，在每个角度区域内对间断有限元基函数积分即可得到初始离散方向对应的求积权重。随着离散方向数的增加，权重系数可能出现负值，通过调整离散方向位置迭代优化求积组以保证权重非负。最后根据二十面体的旋转与对称特性产生全角度区域内求积组。

图1 间断有限元求积组构造流程Fig.1 Construction of discontinuous finite element quadrature sets

图2 平面三角形到球面三角形的投影Fig.2 Projection from flat triangle to spherical triangle

1.2 初始离散方向与基函数

a——线性间断有限元求积组；b——二次间断有限元求积组图3 各次间断有限元求积组Fig.3 Discontinuous finite element quadrature sets of different orders

二十面体映射方法通常选取子三角形区域的重心作为初始离散方向，如图3所示。间断有限元求积组通过细化三角形提高求积组的阶数，对于1阶间断有限元求积组，线性间断有限元求积组(ICLDFE)每组有4个离散方向，二次间断有限元求积组(ICQDFE)每组有9个离散方向。N阶线性间断有限元求积组和二次间断有限元求积组在单个三角形内分别含有4N和9N个离散方向。离散方向数随求积组的阶数呈指数增长，能构造大量离散方向。

选定初始离散方向后，对各已知离散求积点(μm,ηm,ξm)建立基函数bm，其中am,n为基函数系数，m为初始离散方向编号，n为基函数的自由度，与间断有限元基函数的次数有关。线性间断有限元求积组n=4，基函数表达式为：

bm(Ω)=am,1+am,2μm+am,3ηm+am,4ξm

m=1,2,3,4

(1)

二次间断有限元求积组n=9，基函数表达式为：

bm(Ω)=am,1+am,2μm+am,3ηm+am,4ξm+

(2)

在离散角度区域内，基函数在其对应离散方向上的值为1，在其余方向上值均为0。以二次间断有限元基函数建立矩阵方程(3)，通过求逆运算即可求解间断有限元基函数中的未知系数。

(3)

1.3 求积权重及优化策略

对每个离散方向上的基函数在其相应离散角度区域ΔΩ内积分得到求积权重：

(4)

由于上述积分边界是曲线，数值上难以求解。因此，利用雅各比行列式将球面积分(μ,η,ξ)转变为平面积分(u,v)。根据雅各比行列式定义，得：

(5)

其中，|J|为雅各比行列式，即：

(6)

由于坐标系选择的特殊性：

(7)

式(6)可简化为：

(8)

由二十面体的几何特性可得：

(9)

(10)

其中，h为平面三角形的高，由二十面体的几何性质决定：

(11)

将式(9)和(10)代入式(8)，并化简得：

(12)

实际上，随着求积组阶数的增加，中心三角形的权重逐渐减小，从而导致权重出现负值。根据4个初始离散方向得到的线性基函数如图4所示，可看出基函数在相应离散区域并不完全为正，因此需通过调整离散方向的位置以保证权重非负。如图5所示，规定每个离散方向仅在其与平面三角形中心点的连线上移动，且不超过其所在子三角形。定义每个离散方向的比例系数：

(13)

其中：d1为三角形中心到顶点的距离；d2为三角形中心到离散方向点的距离。通过迭代计算调整比例系数使中心三角形的权重等于其对应球面三角形的面积。

已知一个三角形上的求积点之后，根据二十面体的几何特性通过旋转与对称即可得到其余19个面的求积点及权重。图6为2阶线性间断有限元求积组分布图，共有320个离散方向均匀分布球表面，最大权重与最小权重之比为1.58。另外，二十面体的不同三角形面可采用不同阶数或不同次数的求积组，易于实现局部角度细化。

图4 线性间断有限元基函数图像Fig.4 Linear discontinuous finite element basis function

图5 线性间断有限元求积组优化策略Fig.5 Optimization strategy of linear discontinuous finite element quadrature sets

图6 2阶二十面体线性间断有限元求积组Fig.6 Second-order discontinuous finite element quadrature sets on icosahedron

2 精度分析及IRI-TUB基准验证

2.1 球谐函数积分精度

为保证粒子平衡，要求离散求积组在全角度区域内尽可能准确积分高阶球谐函数。但实际屏蔽问题中，局部角度区域内的角通量密度通常是剧烈变化甚至间断的。为研究二十面体间断有限元求积组在局部角度区域内的球谐函数积分精度，用如下公式计算各求积组在1/20

球面上的多项式数值积分：

(14)

其中，a、b、c为多项式阶数，非负整数。表1列出1/20球面内各求积组数值积分结果与数值解的相对偏差。其中ICLDFE-SN为N阶二十面体线性间断有限元求积组；ICQDFE-SN为N阶二十面体二次间断有限元求积组。图7示出了分别在不同多项式阶数下ICLDFE求积组的收敛情况，其中纵坐标为相对误差的绝对值，定义横坐标为平均角度网格步长h：

(15)

由表1可看出，ICLDFE求积组在多项式(0,0,0)与(1,0,0)处的相对误差均为0.00%，ICLDFE求积组能在1/20球面内准确积分0阶以及1阶的球谐函数。同理，ICQDFE求积组则可准确积分2阶及2阶以下的球谐函数。ICLDFE-S4及ICLDFE-S5对于不同多项式阶数的数值积分相对误差均为0.00%。数据表明该求积组在局部区域内具有较高的数值积分精度，随着求积组阶数的增加，数值积分精度逐渐提高。图7表明，平均角度网格步长每减少约1个量级，5阶及5阶以下多项式数值积分的相对误差均下降约4个量级，ICLDFE求积组具有4阶收敛性。

2.2 IRI-TUB基准题

1) 基准题描述

IRI-TUB实验孔道模型为布达佩斯技术大学功率为100 kW的研究反应堆[12-13]，其目的是研究粒子在孔道内的输运行为及验证粒子输运计算方法。由于孔道的存在会使某些区域的角通量密度各向异性分布趋势加剧，当离散求积组数值积分角通量密度精度较低时会导致射线效应问题[14]，极大地影响了屏蔽计算的精度及辐射屏蔽设计的可靠性。因此，采用二十面体线性间断有限元离散求积组进行计算以验证该求积组的数值积分精度及对于实际屏蔽模型的适应性。

表1 1/20球面内各求积组数值积分相对误差Table 1 Relative error of numerical integration for quadrature sets in one-twentieth sphere

图7 不同多项式阶数下数值积分相对误差随平均角度网格步长的变化Fig.7 Relative error of different order polynomials integration versus angular mesh length

图8 IRI-TUB基准题几何结构Fig.8 Geometry structure of IRI-TUB benchmark problem

几何模型如图8所示。堆芯由燃料组件、水及石墨组成。堆芯活性区由24个7.2 cm×7.2 cm的燃料组件组成，高度为50 cm。堆芯外侧配备大型辐射通道，通道内由混凝土填充，混凝土块内有一内径为11.8 cm的圆柱形空腔孔道，孔道外侧包围着厚度为4.5 mm的不锈钢层，总长度为187 cm。沿孔道中心距离孔道入口分别为0、67、121、148、175 cm处装有探测器，可测得孔道内不同位置中子反应率的实验值。

网格划分为3 cm×3 cm×3 cm，差分方式选择菱形差分置零修正方法，多群截面数据由处理172群的MUSE1.0截面数据库得到，考虑各向异性P3阶勒让德多项式近似，求积组分别选用具有相近方向数的PNTN-S60求积组(3 720个方向)与线性间断有限元局部细化求积组(3 584个方向)，求积组分布如图9所示。边界条件均为真空，内迭代收敛准则5×10-4。采用三维离散纵标屏蔽程序ARES[15]直接对该模型进行三维输运计算。

图9 二十面体线性间断有限元局部细化求积组Fig.9 Discontinuous finite element local refined quadrature sets on icosahedron

2) 结果分析

115In(n,n′)115Inm反应率的实验值与计算值比较列于表2。其中C/E(ICLDFE)为采用二十面体间断有限元局部细化求积组的计算结果与实验测量值的比值，C/E(PNTN)为采用60阶勒让德-切比雪夫求积组的计算结果与实验测量值的比值，C/E(文献)为文献中的计算值与实验测量值的比值。文献中的计算值是采用蒙特卡罗程序MORSE-SGC/S接续二维离散纵标程序DOT3.5计算得到。反应率计算公式为：

(16)

表2 IRI-TUB基准题反应率的实验值与计算值比较Table 2 Comparison of measured and calculated reaction rates for IRI-TUB benchmark problem

由表2可看出，采用三维模型直接进行输运计算得到的结果均优于文献中给出的计算结果。对于PNTN-S60求积组，计算结果与实验测量值的相对偏差在35%以内，且最大偏差位于孔道出口附近。这是由于在孔道内部，中子通量密度的各向异性随距孔道入口距离的增加而增强，PNTN-S60仍不能准确描述孔道内中子通量密度的各向异性分布情况，会严重低估中子通量密度。采用二十面体间断有限元局部细化求积组的计算结果更接近实验测量值，与实验测量值最大相对偏差为25%。且随着距孔道入口距离的增加，C/E由0.75增加到1.04。二十面体线性间断有限元局部细化求积组对沿孔道方向的角度区域进行局部细化，提高局部求积精度。因而孔道内的计算值与实验值偏差较小。

图10 IRI-TUB基准题中子能谱的实验值与计算值比较Fig.10 Comparison of measured and calculated spectra for IRI-TUB benchmark problem

图10示出了测量点1和测量点4处大于10 eV中子能谱的计算值与实验值。其中测量点1为孔道入口处，测量点4为距孔道入口148 cm处。从图10可看出，能量为10 eV以上的中子能谱计算值与实验测量值呈现相同变化规律且符合较好。测量点1和4的能量10 eV以上的中子通量密度测量值与实验值之比分别为0.910和0.867。随着距离孔道入口的增加，中子通量密度逐渐衰减，从测量点1至4下降约两个数量级。

3 结论

基于间断有限元思想，构造二十面体线性和二次间断有限元求积组，并对求积方向及权重进行优化。该求积组具有产生大量离散方向且易于局部细化和权重非负的特点。数值结果表明，在1/20球面内，二十面体求积组能准确积分对应阶数以下的所有阶球谐函数，且具有4阶收敛性；IRI-TUB基准题中孔道内各关键点反应率的计算值与实验值的相对偏差在25%以内，中子能谱的计算值与实验值符合较好。初步验证了二十面体间断有限元求积组对于实际深穿透孔道屏蔽问题的适应性与可靠性。但由于二十面体的几何特性，该求积组难以处理反射边界条件。另外，为保证精度同时提高计算效率，生成局部细化求积组的自适应算法有待进一步研究。