改进泄漏项分割技术的二维/一维耦合方法

2021-02-10孔勃然朱凯杰刘保坤张汉郭炯郝琛侯捷李富

哈尔滨工程大学学报 2021年12期

孔勃然，朱凯杰，刘保坤，张汉，郭炯，郝琛，侯捷，李富

(1.清华大学核能与新能源研究院，北京 100084； 2.哈尔滨工程大学核科学与技术学院，黑龙江哈尔滨 150001)

二维/一维耦合方法是一种高保真的一步法计算方法，它将三维方程分解为径向二维方程和轴向一维方程，并通过泄漏项将二者进行耦合求解。受限于计算速度和内存，直接三维全堆芯一步法计算尚存难度[1]，具有相同精度并节省计算量和内存需求的二维/一维耦合方法目前更具实际的工程价值和研究价值，是近年的研究热点[2]。目前采用二维/一维方法的典型程序有DeCART[3]、MPACT[4]、NECPX[5]，HNET[6]、KuaFu[7]等，这些程序普遍采用特征线方法(method of characteristics，MOC)用于二维计算，一维轴向方程则采用扩散方法如节块法(nodal expansion method，NEM)和SP3方法或者输运方法如离散纵标法，泄漏项则采用各向同性或各向异性泄漏项。数值实验表明，采用轴向SN计算并配合各向异性泄漏项的方法精度更高，但与此同时计算量提高、收敛性变差[8]。为了解决计算量大幅提升的问题，MPACT[9]提出方位角傅里叶级数展开的各向异性泄漏项代替完全各向异性的泄漏项，并证明1阶或者2阶傅里叶级数展开就可达到完全各向异性的精度。由于采用了各向异性泄漏项，方程的收敛性相较于采用各向同性泄漏项的方法变差。在迭代过程中，由于泄漏项的存在，二维输运计算过程中可能会计算出负角通量，并有可能导致方程发散。为了解决这一问题，MPACT提出泄漏项分割技术，但是该方法存在计算误差较大的问题[10]。

为了保证精度的同时降低计算负担，本文采用方位角傅里叶级数展开的各向异性泄漏项和轴向SN的组合的二维/一维计算方法，并提出了一套适用于该组合的改进的泄漏项分割技术。该方法不会增加内存负担的同时，相比于传统方法具有更好的效率和收敛性。

1 二维/一维耦合方法

1.1 二维/一维耦合方法的理论

直角坐标系下，稳态形式离散的三维玻尔兹曼输运方程形式为：

(1)

式中：μ为极角；αm为方位角；Σt,g为总截面；φg,m表示第g群m方向上的角通；Qg,m为总源项，包括各向同性的散射源和裂变源。

二维/一维方法将堆芯沿轴向分为若干层，对方程(1)沿轴向积分取平均即可获得径向二维方程：

(2)

(3)

同理，二维/一维方法将堆芯沿径向分为若干个平源区/栅元，对方程(1)沿径向积分取平均即可获得轴向一维方程：

(4)

(5)

从式(2)和(4)可以看出，二维/一维方程通过泄漏项进行耦合，在计算前泄漏项需要从一维/二维方程中获得。通过反复迭代计算二维/一维方程直至特征值和标通量收敛，其原理见图1。

图1 二维/一维耦合输运计算方法Fig.1 2D/1D coupling transport calculation

1.2 方位角傅里叶级数展开的各向异性泄漏项及轴向一维方程

从式(5)知，若采用各向异性泄漏项，虽然轴向方程在空间维度是一维的，但它仍然与方位角有关，此时轴向计算仍需计算所有角度。若轴向方程采用输运方法求解，其计算量相比于各向同性泄漏项显著增加，以KUCA和OECD C5G7-3D RB基准题为例，轴向计算量约占总计算量的20%，可见此时轴向计算非常耗时。为了解决轴向计算量大幅提升的问题，采用方位角傅里叶级数展开的泄漏项代替完全各向异性泄漏项，以实现用少量阶数的轴向方程代替全角度的轴向方程，将角通量和径向泄漏项采用N阶傅里叶级数展开为：

(6)

(7)

将式(6)和(7)代入方程(4)中计算后，对方位角进行积分即可得到零阶、sin阶、cos阶方程：

(8)

(9)

(10)

至此，通过方位角傅里叶级数展开泄漏项，可以将原有全方位角的SN计算转化成 (2N+1)个SN方程进行计算。

2 泄漏项分割技术

2.1 传统的泄漏项分割技术

由式(2)可知，当泄漏项过大时右端源项为负，可能导致输运计算的角通量为负，进而导致算法不收敛。为了解决这一问题，需要采用泄漏项分割技术进行修正，其原理为当右端源项为负值时，修改方程(2)为：

(11)

其中：

(12)

或者：

(13)

其中：

(14)

2.2 改进的泄漏项分割技术

负源项在二维/一维算法中是很难避免的，从算法的机理而言，即使收敛时负源项也可能存在。一定程度的负源项在计算过程中是可以接受的，而输运计算中负角通量是无意义的，因此新的泄漏项分割技术的触发条件为当输运计算出的角通量为负时进行泄漏项分割。

定义泄漏项：

(15)

(16)

定义极角相同，方位角相反的2个方向为正方向和反方向，此时2个方向的泄漏项为：

(17)

(18)

从式(17)、(18)可以发现，正反2个方向的泄漏项区别只在于奇阶泄漏项的正负，其余0阶泄漏项和偶阶泄漏项均相同，利用此特征设计改进的泄漏项分割技术。

对于正向方程，首先将右端源项中奇阶泄漏项去除，若此时源项大于零，则此时不对方程进行修改；反之采用类似于传统泄漏项修正的方法将右端源项置零，左端引入修正截面。具体操作如下：

1)对于正方向，当MOC计算的出射角通量小于0时，进行分割(忽略位置和角度标识)：

(19)

(20)

2)对于反方向，首先在负方向总源项上加上反向修正项，此时正反2个方向的泄漏项相同，相当于对正反2个方向的泄漏项做各向同性处理。

若此时出射角通量仍然小于0，进行和传统修正方法类似分割：

(21)

(22)

该修正方法充分利用了泄漏项的特性，当正方向出现负角通量时，通过上述处理，相当于对正反2个方向的泄漏项做各向同性处理。虽然在特定方向上的新方程不等价于原方程，但对所有方向新修正的方程进行积分，新方程积分等价于原方程积分，这就保证了二维计算和一维计算得到的标通量一致等价，因此该修正方法改善了精度和收敛性。

3 精度和收敛性分析

3.1 计算精度分析

1 强泄漏问题

该算例的几何模型如图2所示，径向上燃料棒3×3布置，最中心为慢化剂，径向上采用3个反射边界条件和1个真空边界条件，轴向上采用反射边界条件。轴向上布置2层燃料层和1层慢化剂层。可以将此算例类推为径向7×7、17×17算例，以研究不同泄漏强度对精度的影响。该算例中以完全各向异性泄漏项方法的结果作为基准解，分别测试不采用修正，传统泄漏项修正和新的泄漏项修正方法，其结果如表1所示。数值实验证明采用改进的泄漏项修正方法可以降低修正误差，其结果与不修正方法近似且解决了收敛性的问题。采用修正方法、传统修正、改进修正方法时，三者的迭代次数完全相同，而单次迭代时间也几乎相同，因此此算例中，改进泄漏项方法和传统修正方法的计算时间几乎一致。

图2 强泄漏算例几何模型Fig.2 Geometry of high leakage case

表1 强泄漏算例计算结果Table 1 Calculation results of high leakage cases

2 三维C5G7基准题

C5G7基准题是OECD/NEA发布的带控制棒的非均匀输运计算基准题[11]，分为无棒、半插棒和全插棒。参数选取为全角度64个方位角，8个极角，特征线密度0.01 cm。径向划分为51个栅元，每个栅元划分为5个圆环16个扇区。轴向分为18层，每层3.57 cm。表2为G5G7计算结果，结果显示采用方位角傅里叶级数展开的泄漏项并采用改进的泄漏项分割技术的结果可以达到较好的结果，1阶或者2阶的结果就可以和各向异性的结果相当。图3为C5G73D-Unrodded功率误差，可见采用傅里叶级数展开的泄漏项和各向异性泄漏项的结果几乎一致。

图3 C5G7-Unrodded功率误差Fig.3 Power error of C5G7-Unrodded cases

表2 C5G7计算结果Table 2 Results of C5G7 cases

3.2 收敛性分析

为了分析采用改进的泄漏项修正方法的2D/1D耦合算法的收敛性，选用KUCA基准题[12]进行研究。KUCA是日本东京大学模拟小型反应堆的临界装置，其堆芯包括燃料、慢化剂、控制棒。选取其中KUCA空隙问题(情况1)作为算例进行验证，基准解为0.977 8，其结果如表3所示。从计算结果中可以看出，不采用修正方法时，1阶和3阶傅里叶级数展开方法是发散的，其原因为负角通量过多，导致了负通量的产生，进而导致算法发散。而传统修正方法1阶傅里叶级数展开时仍然发散，其原因为传统方法的所有角度的修正方程积分不等价于原方程，当负角通量过多时，修正后的方程和原方程偏差过大，二维计算统计出的标通量不等价于一维计算统计出的标通量，进而导致整体算法无法收敛。结果证明改进的泄漏项分割方法可以解决负角通量造成的发散问题。表4展示了不同修正方法的迭代次数对比，由于每次迭代计算时间几乎相同，因此迭代次数即反映了计算时间的对比,从结果中可以看出几种方法的计算时间几乎相同。