基于三维特征线方法的矩形CMFD研究
Rectangle CMFD Based on the 3D Method of Characteristics
摘要: 三维特征线方法(the Method of Characteristics, MOC)作为确定论方法,能够满足中子输运求解的计算要求,并且有着优秀的几何适应性、较高的计算精度、简洁易实现的计算模型及强大的并行潜力,但是在高保真、精细化建模的需求下,它对计算量有着巨大的需求,因此需要一个加速方法来提高程序的计算效率,而粗网有限差分方法(coarse mesh finite difference method, CMFD)作为使用的较多的加速方法,能够取得非常显著的加速效果。本次研究基于三维特征线程序ANT-MOC,并使用自建计算模型对CMFD的加速效率和加速精度进行验证,最终结果表明本次研究的矩形CMFD方法在计算组件级的几何模型时能够节省10倍以上的计算时间,且在均匀矩形网格划分的情况下能实现无损加速。
Abstract: As a deterministic neutron transport solving method, the 3D method of characteristics (MOC) has excellent geometric adaptability, high calculation accuracy, simple and easy-to-implement computational model and strong parallel potential. However, under the demand of high-fidelity and refined modeling, it has a huge demand for computation, so it needs an acceleration method to improve the computational efficiency of the program, and the Coarse Mesh Finite Difference method (CMFD) is a widely used acceleration method that can achieve a very significant acceleration effect. In this study, we use ANT-MOC, a 3D MOC program, to count, and self-built calculation model is used to verify the acceleration efficiency and acceleration accuracy of CMFD. The final results show that the rectangular CMFD method can save more than 10 times the calculation time when calculating the component level geometry, and can achieve lossless acceleration in the case of uniform rectangular meshing.
文章引用:杨缙卿, 胡赟, 单浩栋. 基于三维特征线方法的矩形CMFD研究[J]. 核科学与技术, 2025, 13(3): 165-173. https://doi.org/10.12677/nst.2025.133017

1. 前言

ANT-MOC作为一个运用MOC方法求解中子输运方程的程序[1] [2],具有优秀的几何适应性和较高的计算精度,然而在追求较高的计算精度时,特征线程序必须采用足够密集的特征线且对于平源近似区的划分也要足够的精细,因此对计算量有着大量的需求。且由于ANT-MOC程序采用直接三维特征线法来进行中子输运的计算,与使用二维/一维耦合求解三维特征线中子输运方法的程序相比,对计算资源的需求更大,因此迫切需要一个加速方法来加速三维特征线程序ANT-MOC的输运计算。

在诸多特征线程序的加速方法中,粗网格有限差分方法(CMFD)有着巨大的优势。首先,由于CMFD计算和MOC计算收敛的一致性,CMFD方法在理论上能够实现无损加速;其次,就加速效率而言,CMFD方法的加速效果往往能够达到五倍以上,相对于其他加速方法计算效率的提示非常明显[3];最后,CMFD方法在国内外已经进行了大量的研究,计算手段相对成熟,它的差分格式最早由K.S. Smith提出[4],并在后续的几十年内被广泛应用到计算程序之中。

CMFD方法的核心思想是在MOC细网格的基础上构建粗网格并计算粗网格上的等效均匀化参数,然后在此基础上构造一个“等效”的中子扩散问题,通过扩散问题的求解结果来修正MOC细网格上的初始通量,以达到减少迭代次数和计算时间的效果[5]

本次研究使用自建模型对研究的矩形CMFD程序进行加速效果和精度验证,研究结果证明了CMFD方法的加速效果非常显著且计算精度也能达到计算要求。

2. 理论模型

2.1. 特征线方法中的中子平衡关系

中子输运方程的表达形式如下[6]

1 v Φ( r ,E, Ω ,t ) t + Ω Φ( r ,E, Ω ,t )+ Σ t ( r ,E,t )Φ( r ,E, Ω ,t ) = Q S ( r ,E, Ω ,t )+ Q f ( r ,E, Ω ,t ) (1)

上式中,v为中子速度;t表示时刻; r 表示中子位置; Ω 表示中子运动方向角度向量;E表示中子能量; Σ t 表示中子总截面; Φ 表示单位立体角内的中子运行轨迹之和; Q S 表示散射中子源项; Q f 表示裂变中子源项。

在中子输运方程达到平衡时,可以忽略时间对方程的影响。后续对能量进行离散化处理,即划分G个能群,对于能量位于能群g中的中子输运方程可以表示为:

Ω Φ g ( r , Ω )+ Σ t ( r ,E ) Φ g ( r , Ω )= Q S,g ( r , Ω )+ Q f,g ( r , Ω ) (2)

再对空间进行离散化处理,对空间进行分区处理( i=1,2, ),且假设裂变中子的发射方向是各项同性的,那么对空间网格ig能群、角度 Ω 的中子输运方程可以表示为:

Ω Φ g,i + Σ t Φ g,i = Q S,g,i + Q f,g,i = S g,i (3)

在(3)式中, Q S,g,i 代表单元网格内流入流出的中子变化率, Q f,g,i 代表单元网格内各种反应消耗的中子, S g,i 代表单元网格内通过散射和裂变产生的中子。即满足中子平衡方程:

中子变化率 + 反应消耗率 = 源项产生率

2.2. CMFD方法概述

CMFD在加速三维特征线程序过程中,主要分为以下三个步骤,如图1所示:

Figure 1. Schematic diagram of accelerating MOC using CMFD

1. 利用CMFD加速MOC的示意图

1. 压缩过程

通过等效的均匀化将MOC细网格中的截面、扩散系数等归并到粗网格中,可以对中子输运方程在 4π 角度积分:

J g ( r )+ Σ t g ( r ) Φ g ( r )= χ g ( r ) k g =1 G v g ( r ) Σ f g Φ g ( r ) + g =1 G Σ s g g Φ g ( r ) (4)

式(4)中, r 表示空间位置;g表示所处能群示 J g ( r ) 为中子流密度; Σ t g ( r ) 为总截面; Φ g ( r ) 为中子通量; χ 为中子发射光谱, χ g ( r ) 为中子裂变发射概率; v 是每次裂变释放的平均中子数; Σ f g 为中子裂变截面; Σ s g g 表示散射截面。

在任意空间 r V 进行积分,得到式(5):

dr J g ( r )+ dr Σ t g ( r ) Φ g ( r ) = dr ( χ g ( r ) k g =1 G v g ( r ) Σ f g Φ g ( r ) + g =1 G Σ s g g Φ g ( r ) ) (5)

指定CMFD网格为 V C j ,将该粗网格中的细网格指定为 V i ,则有:

dr J g ( r )+ ij dr Σ t i,g ( r ) Φ g ( r ) = ij dr ( χ g ( r ) k g =1 G v g ( r ) Σ f i, g Φ g ( r ) + g =1 G Σ s i, g g Φ g ( r ) ) (6)

在细网格内,可以近似平源区处理,即MOC细网格内材料视为均匀分布。因此,对于任一细网格,认定其体积、截面和中子发射谱等都为常数:

dr J g ( r )+ ij Σ t i,g ( r ) ϕ ¯ i,g V i = ij ( χ i,g ( r ) k g =1 G v i, g ( r ) Σ f i, g ϕ ¯ i,g V i + g =1 G Σ s i, g g ϕ ¯ i,g V i ) (7)

根据高斯定理将中子流项写为曲面积分的形式:

ds J g ( r ) n + ij Σ t i,g ( r ) ϕ ¯ i,g V i = ij ( χ i,g ( r ) k g =1 G v i, g ( r ) Σ f i, g ϕ ¯ i,g V i + g =1 G Σ s i, g g ϕ ¯ i,g V i ) (8)

至此,可以将粗网格中的均匀化常数项利用细网格中的常数表达出:

χ C i,g = ij ( χ i,g g =1 G v i, g Σ f i, g ϕ ¯ i, g V i ) ij ( g =1 G v i, g Σ f i, g ϕ ¯ i, g V i ) (9)

v C g Σ C,f g = ij ( v i,g Σ f i,g ϕ ¯ i,g V i ) ij ( ϕ ¯ i,g V i ) (10)

Σ C,f g = ij ( Σ f i,g ϕ ¯ i,g V i ) ij ( ϕ ¯ i,g V i ) (11)

至此,除了中子流项没有得到处理之外,粗网中其他项都可以用细网格内的已知项通过等效均匀化表达出来。

从整个压缩过程中我们可以看到,CMFD方程的解依赖于通过MOC计算所求得的标通量 ϕ ¯ C j,e ,且CMFD网格上的求解最终也会反馈到MOC细网格输运计算中去,因此在整个CMFD加速计算过程中,理论上在迭代收敛的时候不会产生近似带来的误差,故而两者的解有着一致性。

2. CMFD网格上求解

可通过将上述过程中求解出的等效均匀化参数,代入到CMFD求解方程之中,计算出粗网格上的中子流和通量分布。

首先,将粗网格j的边界面 S j 分解为若干个子界面(假设分为H个),在用 J ˜ j,h,g 表示交界面h流出的中子流时,有以下中子平衡方程:

1 V C j h=1 H J ˜ j,h,g + Σ C,t j,g ϕ C j,g = χ C j,g k g =1 G v C j, g Σ f i, g ϕ C j, g + g =1 G Σ C,s i, g g ϕ C j, g (12)

利用扩散近似进行求解,则通过菲克定律,有:

J ˜ j,h,g A j,h =u( j,h ) D ^ j,h,g ( ϕ C I( j,h ),g ϕ C j,g ) (13)

其中, D ^ j,h,g 是扩散系数, u( j,h ) 表示中子流方向。对于 D ^ j,h,g 的值的求解,假设 Δ r j 为CMFD网格j的中心到边界面h的距离,那么其有着如下的表达形式:

D ^ j,h,g = 2 D j,g D I( j,h ) D j,g Δ r I( j,h ) + D I( j,h ) Δ r j,g (14)

其中, D j,g D I( j,h ) 分别为空间网格j及其边界面h所对应的空间网格内按其内的MOC细网格中子通量平均后的扩散常数。以 D j,g 为例,其表达形式为:

D j,g = ij 1 3 Σ t i,g ϕ ¯ i,g V ij ( ϕ ¯ i,g V i ) (15)

由于扩散方程较输运方程精度偏低,为了提高CMFD的计算精度,做出如下修正:

J ˜ j,h,g A j,h =u( j,h ) D ^ j,h,g ( ϕ C I( j,h ),g ϕ C j,g ) D ˜ j,h,g ( ϕ C I( j,h ),g + ϕ C j,g ) (16)

其中, D ˜ j,h,g 作为修正项中的一部分,它可以以上一次迭代的过程中的数值进行计算:

D ˜ j,h,g = u( j,h ) D ^ j,h,g ( ϕ C I( j,h ),g ϕ C j,g ) J ˜ j,h,g A j,h ϕ C I( j,h ),g + ϕ C j,g (17)

因此,最后的CMFD方程如下:

1 V C j h=1 H A j,h ( u( j,h ) D ^ j,h,g ( ϕ C I( j,h ),g ϕ C j,g ) D ˜ j,h,g ( ϕ C I( j,h ),g + ϕ C j,g ) ) + Σ C,t i,g ϕ C j,g = χ C j,g k g =1 G v C j, g Σ f i, g ϕ C j, g + g =1 G Σ C,s i, g g ϕ C j, g (18)

为了更进一步方法求解,我们对式(18)做进一步的整理,将相同区域内的粗网格标通量整理为一项,则有:

[ V C j Σ C,t i,e + h=1 H A j,h ( u( j,h ) D ^ j,h,e D ˜ j,h,e ) ϕ C j,e h=1 H A j,h ( u( j,h ) D ^ j,h,g + D ˜ j,h,g ) ] ϕ C I( j,h ),g g =1 G Σ C,s i, g g ϕ C j, g = χ C j,g k g =1 G v C j, g Σ f i, g ϕ C j, g (19)

将其写为矩形形式为:

A ϕ C S ϕ C = 1 k M ϕ C (20)

通过对式(20)的求解,便可以求出CMFD粗网格上的平均通量,并且求出有效倍增因子keff

3. 映射并修正

在CMFD方程求解完成后,需要利用CMFD求解得到的标通量分布来更新在MOC方程中求得的细网格通量。一般采取前后粗网的比值作为更新系数:

ϕ i new = ϕ C,j new ϕ C,j old ϕ i old (21)

在(21)式中, ϕ i old 指MOC细网通量更新前的通量, ϕ C,j old 指从MOC计算得到的结果 ϕ i old 中通过空间平均和能群求和得到的CMFD粗网格通量, ϕ C,j new 指通过迭代后求得的CMFD粗网格通量。

将粗网格的计算结果映射到细网格中,用于修正MOC细网格的参数,为MOC中子输运方程提供良好的初值,从而减少迭代次数。

综上,CMFD方法加速三维特征线程序的手段主要是通过加速MOC的内迭代(源迭代)的收敛加速中子输运的计算过程,主要通过两种机制来加速MOC输运计算,首先是减少迭代次数,通过粗网格的快速计算,为MOC细网格提供更接近真实解的初始值,从而减少MOC计算的迭代次数;最后,是简化计算降低计算规模,在方法上使用计算起来更为简化的扩散近似来代替输运方程,以此降低CMFD网格上的计算量,在规模上使用更为粗糙的CMFD网格,通过等效均匀化将MOC网格的材料属性以及通量压缩到CMFD网格中。因此在整个计算过程中,求解CMFD方程要比直接进行MOC输运求解的效率高很多,且在CMFD加速方法中,也可以通过对能群进行压缩来更进一步的提高MOC程序的计算效率。

另外,使用CMFD加速方法的迭代特性使得该方法在计算精度上不会降低,且能够实现无损加速,即每次CMFD迭代的初值都会依赖MOC的标通量,且CMFD上每次迭代计算之后的数据都会以映射的方式反馈到MOC的输运计算之中,因此尽管CMFD方法求解的是扩散方法,但是计算精度仍然达到了输运方程的精度要求。

2.3. 计算流程

图2是CMFD方法在MOC程序中的应用流程图,可以看到在CMFD加速MOC内迭代之后,会将所计算得到的结果反馈到MOC细网格上,并且会再次在MOC细网格的基础上进行一次内迭代的求解,因此使用CMFD方法加速MOC输运计算的计算精度并不会因为CMFD方法是求解的中子扩散方程而下降。相反,最终因为会再次回归MOC输运计算来进行求解,因此使用CMFD方法加速MOC输运计算在精度上是无损的。

3. 实验验证

自建栅元结构如图3所示,其中不同材料以颜色进行区分,从内至外依次为燃料、空隙、包壳、冷却剂。各尺寸已进行标注,燃料半径0.4096 cm,包壳内半径0.4178 cm,包壳内半径0.4750 cm,整体栅元为边长为1.26 cm的正方形,轴向高度2 cm。

组件示意图如下。

自建堆芯组件结构如图4所示,以栅元正方形结构为基础进行周期性排列,共计289个栅元结构,增加的组件壁结构材料几何尺寸如图所示,组件内壁对边距为21.42 cm,组件外壁对边距为21.66 cm,轴向高度为2。

Figure 2. CMFD flow chart

2. CMFD流程图

Figure 3. Schematic diagram of Grid

3. 栅元示意图

Figure 4. Schematic diagram of component

4. 组件示意图

以RMC计算值作为参考值,则对于单栅元和单组件,RMC的计算值如下表1

Table 1. Calculation of RMC

1. RMC计算

计算对象

Keff

标准差

单栅元

7.7551E−02

6.6E−06

单组件

1.421650

3.39E−05

基于该模型,使用特征线程序ANT-MOC及其CMFD模块进行加速效果验证,结合模型的实际情况,采用64方位角进行计算,特征线间隔划分为0.01cm使用6极角,间距设置0.75 cm。

首先记录该特征线程序在不使用CMFD加速的情况下,其keff、迭代次数和计算耗时如下:

Table 2. Calculation without CMFD

2. 不开启CMFD计算

计算对象

keff

迭代次数

计算时间/S

单栅元

0.076233

173

1.73E+02

单组件

1.41894

541

2.12E+04

可以看到使用纯粹的ANT-MOC输运计算,计算结果与基准值之间的误差在可接受范围内。开启CMFD模块进行计算,其计算结果如下:

Table 3. Calculation with CMFD

3. 开启CMFD计算

计算对象

keff

迭代次数

计算时间/S

时间加速比

单栅元

0.076212

31

2.78E+02

0.622

单组件

1.41972

26

1.62E+03

13.08

表2表3的计算结果中可以看出,不论是栅元的计算还是组件的计算,CMFD加速方法都能够明显的减少计算时的迭代次数,并且计算规模越大,迭代次数减少的越明显。

另外,随着计算对象规模的增大,CMFD方法起到的加速效果越明显,在计算规模过小时,CMFD方法不仅不会起到加速效果,甚至会增加计算时间,这是因为CMFD方法本身的压缩映射过程和扩散计算也需要耗费一定的时间。

CMFD方法的最终计算结果与MOC输运计算结果的偏差较小,在本次实验中,在栅元量级的计算偏差为28 pcm,在组件量级的计算中,偏差为55 pcm,这与前面提到的CMFD能够实现无损加速相吻合。

4. 结论

本次研究建立了矩形CMFD方法的理论模型,并通过自建模型对均匀划分的矩形CMFD方法进行了效果验证,其结果如下:

对于均匀矩形网格划分的CMFD,在自建模型的栅元和组件计算中,迭代次数都明显减少,其中计算规模越大,迭代次数减少的越多。

就计算时间来说,栅元量级在使用CMFD加速后计算时间更长,这是由于迭代过程中CMFD模块本身就要占用一定的时间,因此对于过小的几何对象CMFD的加速效果非常差,本次实验出现的现象与理论相符;而对于组件量级的计算时间得到明显减少,加速效果达到了10倍以上,并实现了无损加速。综上可以证明本次研究的CMFD在均匀矩形网格划分的情况下取得了巨大的成功。

NOTES

*通讯作者。

参考文献

[1] 单浩栋. 面向数值反应堆通用几何的三维特征线输运算法研究与并行程序开发[D]: [博士学位论文]. 北京: 中国原子能科学研究院, 2020.
[2] 单浩栋, 徐李, 胡赟, 胡长军, 王根, 方雅. 三维高保真特征线程序ANT-MOC的开发与测试[J]. 核动力工程, 2020, 41(z1): 1-5.
[3] 高明敏, 曹欣荣, 朱成林, 王江萌. 二维MOC的CMFD加速计算模块开发与验证[J]. 原子能科学技术, 2017, 51(1): 106-112.
[4] Gunow, G.A. (2018) Full Core 3D Neutron Transport Simulation Using the Method of Characteristics with Linear Sources. Massachusetts Institute of Technology.
[5] Zhong, Z., Downar, T.J., Xu, Y., DeHart, M.D. and Clarno, K.T. (2008) Implementation of Two-Level Coarse-Mesh Finite Difference Acceleration in an Arbitrary Geometry, Two-Dimensional Discrete Ordinates Transport Method. Nuclear Science and Engineering, 158, 289-298.
https://doi.org/10.13182/NSE06-24TN
[6] 谢仲生, 尹邦华, 潘国品. 核反应堆物理分析下册[M]. 北京: 原子能出版社, 1994.