二维中子输运问题的离散纵标短特征线方法研究

刘 聪，张 斌，张 亮，包博宇，赵好强，陈义学

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

离散纵标法(SN)是粒子输运数值计算常用的确定论方法之一[1]，广泛应用于堆芯物理分析和核设施屏蔽计算。空间变量离散是SN方程数值求解的重要组成部分，常见的离散方法包括有限差分法、有限元法等，不同离散方法各有优劣。空间离散方法的正定性、截断误差、网格步长敏感性等影响着输运计算的可靠性。有限差分方法具有成熟的数学基础且是各类SN程序普遍应用的空间离散方式，菱形差分(DD)[2]是最具代表性的低阶差分方法，DD在一维平板几何中具有二阶收敛精度，但有负通量问题且在粗网格下计算精度较低。由DD发展出多种非线性低阶差分方法，包括负通量置零修正菱形差分(DZ)、带权重差分(WD)[3]、定向θ权重差分(DTW)[4]、指数定向差分(EDW)[5]和指数迭代差分(EDI)[6]等，这些差分方法通过引入修正或非线性假设保证角通量密度非负，在一定程度上改善了DD非物理振荡和粗网精度低的问题。有限元法[7-8]利用基函数在网格内将角通量密度进行展开，求解离散后的弱形式输运方程，具有较高计算精度和较好数值扩散特性，但多维情况下随着展开阶数升高，负通量抑制技术[9]并不能严格保证全局角通量密度非负。Lathrop[10]提出SN方程的常数短特征线离散方法，假设网格内源强分布构造半解析的特征线解，通过输运方程的空间矩守恒推导网格角通量密度分布。Dehart等[11]将其推广至二维任意形状网格，开发了NEWT输运计算程序。Larsen和Alcouffe[12]以及Mathews等[13]分别基于结构和非结构网格研究了线性短特征线方法，后者对负通量不进行修正。按照对网格源强分布处理方式的不同，自1990年至今陆续有多种基于短特征线离散的数值方法被提出。

相比差分方法仅考虑中子角通量密度沿单个空间变量的变化，短特征线方法基于半解析的中子衰减解，根据离散角度考虑网格出射边界与入射的贡献关系，减弱了因空间离散误差造成的非物理振荡。常数短特征线格式天然非负但计算精度不足，本文重点研究SN方程的线性短特征线空间离散方法，改善粗网精度、提高计算效率，通过对角通量密度分布的线性斜率进行旋转修正，严格保证源强和通量密度非负，避免非物理结果。

1 短特征线空间离散方法

1.1 SN短特征线方法与长特征线方法的比较

图1 SN短特征线和MOC方法的扫描示意图Fig.1 Sketch of mesh sweep for SN short characteristic and MOC methods

短特征线方法是SN方程的空间变量离散方法，区别于组件、堆芯分析中常用的特征线方法(MOC方法，也称为长特征线方法)。SN短特征线和MOC方法均是基于微分-积分形式中子输运方程的确定论求解方法，方程中连续角度变量直接离散为1组有限离散方向，两者最主要的理论差异体现在空间变量离散和网格扫描求解上。SN短特征线方法求解基于中子守恒的网格平衡方程，根据空间离散格式的理论假设，由网格入射角通量密度、网格源强求解得到网格平均和出射角通量密度，按离散方向进行角度-网格输运扫描。短特征线离散根据特征线解与SN方程的空间矩守恒，推导网格入射与出射角通量密度的数值关系。MOC方法由平源或线性源近似对输运方程进行积分可得到特征线解。每个离散方向下划分1组平行射线穿过计算系统，几何运算求出射线与网格交点、线段长度，利用射线追踪技术进行扫描计算，求解各条射线在网格内入射、出射和该段平均角通量密度。SN短特征线和MOC方法按图1所示路径对计算模型进行扫描。短特征线方法基于SN网格平衡方程，针对局部网格，而MOC方法是沿穿过整个计算系统的特征线进行求解，针对网格内各段子射线。此外，SN短特征线和MOC方法在网格平均通量密度计算、边界处理等方面存在差异。由于计算问题的不同，对角度求积组、各向异性散射展开阶数的选择也不同。

1.2 常数短特征线空间离散方法

稳态中子输运方程中的角度变量采用SN方法直接离散，能量变量分群处理，二维xy几何下的SN方程如下。

Σt,gψm,g(x,y)=Qm,g(x,y)

(1)

式中：ψ为中子角通量密度；Σt为宏观总截面；x和y为空间变量；Q为源强，包含散射源、裂变源及外源；μ和η为角度向量与坐标轴夹角余弦；m和g分别为离散方向编号和能群编号，后文为简洁表示，省略离散方向和能群编号下标。

根据离散方向导数关系将式(1)改写为射线方程(式(2))，该方程具有式(3)形式的解。

(2)

(3)

常数短特征线方法也称步特征线(SC)方法，假设网格源强和边界入射角通量密度为常数分布。源强由初值或前次迭代结果给出，入射角通量密度由边界条件或上风向网格出射值给出。建立计算网格的局部坐标系，二维几何下角特征线与网格的相交情况如图2所示，角特征线上、下分别由左边界和下边界入射贡献。

a——交于上边界；b——交于右边界图2 角特征线与网格关系图Fig.2 Relationship between angular characteristic line and mesh

(4)

通过特征线解与目标函数的矩守恒可推导SC方法的边界出射值与网格平均值，其中目标函数为常数分布，仅需零阶矩守恒(式(5))。

(5)

(6)

(7)

(8)

计算网格的出射值作为下风向网格入射值，完成角度-网格扫描，更新网格通量矩和散射源，采用源迭代方法进行求解。

1.3 线性短特征线空间离散方法

SC方法具有正定性，在源强和入射值非负的条件下可保证角通量密度计算值非负，但由于基于常数假设，在粗网格下计算精度较低。线性短特征线(LC)方法基于线性假设改善了粗网格的计算精度，降低了输运计算对网格步长的敏感性。计算网格局部坐标系下假设网格源强、入射角通量密度为线性分布，式(2)的SN方程改写为式(9)形式。

0

(9)

(10)

Ψabove(x,y)=

Ψbelow(x,y)=

需要基于特征线解构造线性分布的出射ψi,out、ψj,out和网格ψcell角通量密度分布函数。

(11)

式中，ψx和ψy为角通量密度在网格内分布函数的线性斜率。

(12)

图3 网格角通量密度求解示意图Fig.3 Sketch for solving mesh angular flux

基于零阶矩守恒可求解式(11)和(12)中边界角通量密度平均值。

(13)

(14)

(15)

(16)

目标函数为线性分布，边界线性斜率按式(17)根据差商近似求解。

(17)

对于相交边界分别计算各段斜率。

(18)

根据一阶矩守恒，将式(12)代入构造边界的整体线性分布函数。

(19)

最终得到相交边界的角通量密度分布函数的线性斜率。

(20)

(21)

(22)

(23)

按标准球谐函数展开法计算散射源，更新网格源强进行迭代计算。LC方法不是天然非负的，当网格步长过大时角通量密度分布可能在边界或网格内出现负值，由于负通量是非物理的同时可能导致网格负散射源，因此进行旋转斜率修正强制边界或网格线性分布最小值为零，该修正方法破坏了网格一阶矩守恒关系，但避免了数值求解的非物理结果，减弱了负通量对迭代的扰动。短特征线离散格式SC和LC已嵌入多维SN输运计算程序ARES[14]中的二维求解模块DONTRAN2D。

2 数值验证

2.1 Gelbard固定源问题

该测试例题是双群均匀介质固定源问题，取自美国阿贡国家实验室基准题手册ANL-7416[15]，四周为反射边界，文献给出P19阶球谐函数解作为参考基准。模型几何如图4所示，材料截面和源强信息列于表1，表中，ν为每次裂变释放的中子数，Σf为裂变截面，Σs为散射截面。

图4 Gelbard固定源问题几何示意图Fig.4 Geometric model of Gelbard fixed-source problem

gνΣf,g/cm-1Σs,g-1/cm-1Σs,g-2/cm-1Σt,g/cm-1源强/(cm-3·s-1)10.06.947×10-32.343 4×10-29.210 4×10-26.546 0×10-320.00.04.850×10-31.008 77×10-11.770 1×10-2

网格尺寸设为1 cm，角度求积组选取S16阶EON求积组[16]，通量密度收敛准则为5×10-5。沿y=80和y=139的计算结果示于图5，SC和LC的计算结果与参考解吻合良好。强吸收介质引起SN方法的射线效应，造成沿y=139处的通量密度结果一定程度的振荡。

图5 各能群中子通量密度计算结果Fig.5 Neutron flux results for each group

2.2 自设固定源问题

图6 单群固定源问题几何模型Fig.6 Geometric model of one-group fixed-source problem

为分析各空间离散方法对屏蔽问题的计算精度，自设单群均匀介质固定源问题。模型如图6所示，左边界和下边界为反射边界，其余为真空边界，材料总截面为0.5 cm-1，散射比为0.5，该模型长度为15倍中子自由程。使用多群蒙特卡罗程序MCMG[17]对该问题进行模拟，沿下边界设置6个计数器，每个边长为5 cm，统计区域平均中子通量密度作为参考解。

图7 单群固定源问题的各区域平均中子通量密度和相对偏差Fig.7 Regional average neutron flux and relative deviation for one-group fixed-source problem

图7示出了在网格步长为0.2 cm条件下SC和LC的计算结果以及MCMG各计数器的统计误差。可看出LC结果与MCMG基准的相对偏差小于1.4%，SC计算精度不及LC方法。

网格步长分别设为0.1、0.2、0.5、1、2.5和5 cm，对短特征线、差分空间离散方法进行测试，由各区域相对偏差按式(24)计算离散方法的均方根(RMS)偏差。

(24)

式中，φn,calc和φn,ref分别为第n个区域的平均中子通量密度计算值和基准值。

图8示出了各方法均方根偏差关于网格光学厚度(ΣtΔl，Σt为材料总截面，Δl为网格边长)和计算时间的变化情况。可见，LC和EDW对于该问题的计算效率最高，在相近时间内可获得更高的计算精度，对网格步长敏感性较低，网格步长为2.5倍自由程时的RMS偏差约为10%。但EDW在网格步长为0.05和0.1倍自由程时计算未收敛，由于该方法假设通量密度在网格内为指数函数分布，网格划分过细时网格内的通量密度衰减很小，导致指数因子拟合出现问题。SC方法基于常数假设，计算结果对于网格步长十分敏感，粗网下精度不足。

图8 各空间离散方法的RMS偏差分布Fig.8 Distribution of RMS deviation for each spatial discretization method

2.3 Stepanek沸水堆栅元临界基准题

Stepanek等[18]提出了轻水堆基准题，用于比较不同确定论方法的计算精度与效率，本文选取其中沸水堆(BWR)栅元问题进行计算分析。该模型为双群临界问题，内部均匀化燃料棒由轻水包围，四周为反射边界，几何模型如图9所示，材料截面列于表2，裂变谱快群为1.0、热群为0.0。

图9 沸水堆栅元几何示意图Fig.9 Geometric model of BWR cell problem

材料gνΣf,g/cm-1Σs,g-1/cm-1Σs,g-2/cm-1Σt,g/cm-1116.203×10-31.780×10-11.002×10-21.966 47×10-121.101×10-11.089×10-35.255×10-15.961 59×10-1210.01.995×10-12.188×10-22.220 64×10-120.01.558×10-38.783×10-18.878 74×10-1

网格步长约0.2 cm，特征值收敛准则设为10-6，参考解来自MCNP的蒙特卡罗计算结果[19]。特征值Kinf和区域平均通量密度的计算结果列于表3，同时给出了ARES的有限差分离散的计算结果，区域平均通量密度按燃料区快群结果进行了归一。该网格划分下几种空间离散格式的特征值结果与参考值的偏差均小于100 pcm。其中LC结果与参考值最接近，仅偏大16 pcm，其次是DZ和WD格式较准确。

改变网格步长进行网格敏感性的测试，随网格数量和计算时间改变的特征值结果变化示于图10。可见LC在各网格划分下的结果均优于其他离散格式，6×6网格划分时偏差也仅为58 pcm。几种差分离散方法对于该问题的计算效率差异较小，DZ与WD基本一致且略优于EDW和DTW。SC随着网格步长的增大计算结果出现较大偏差，计算精度与网格步长存在较强关联。

表3 沸水堆栅元问题的计算结果Table 3 Result of BWR cell problem

图10 沸水堆栅元问题特征值Fig.10 Eigenvalue result of BWR cell problem

3 结论

本文采用离散纵标短特征线方法求解二维中子输运方程，研究了基于常数和线性分布的短特征线空间离散方法，并将该方法应用于多维SN输运计算程序ARES中的二维求解模块DONTRAN2D中。固定源和临界测试例题的结果表明，线性短特征线方法对网格步长敏感性较低，相比常数短特征线和低阶差分离散具有更高的粗网精度。未来将扩展至三维几何并在实际复杂工程问题上进行应用。