三维MOC的球床堆芯CMFD加速研究

文宇晨， 郝琛,2， 郭建， 王毅箴

(1.哈尔滨工程大学 核安全与仿真技术重点学科实验室，黑龙江 哈尔滨 150001；2.中国核电工程有限公司， 北京 100840；3.清华大学 核能与新能源技术研究院， 北京 100084)

目前球床式反应堆的物理分析软件有：基于中子扩散理论的VSOP[1]；在VSOP基础上优化后采用更先进两步法进行高温气冷反应堆物理和燃料循环计算的PANGU[2]；美国首选球床堆物理分析软件PEBBED[3]能够在中子扩散计算的迭代过程中求解稳态情况下的中子通量、燃耗、功率等反应堆参数。但对反应堆堆芯提供精细的通量分布可以降低保守安全裕量估计，提高反应堆经济效益，因此对球床式反应堆进行精细化中子输运计算是十分必要的。

郭建等[4]针对球床反应堆特殊的堆芯结构，以构造实体几何方法(constructive solid geometry, CSG)为基础，建立精细的球床反应堆模型，并用三维特征线方法(three dimension-method of characteristics, 3D-MOC)求解中子输运方程，开发了准确描述球床高温气冷堆几何的三维特征线程序(method of characteristics for pebble-bed HTR, MOCP)。为解决3D-MOC方法计算时间长和迭代次数多等问题，MOCP使用稀疏长特征线方法(macro-track transport acceleration method, MTTA)作为3D-MOC的加速算法。MTTA相较基于中子扩散方程的加速算法需要更大的计算量，而基于中子扩散方程的粗网有限差分方法[5](coarse mesh finite difference, CMFD)加速效果好且便于实施，通常作为轻水堆精细化输运计算的加速算法。但传统CMFD可能出现耦合扩散系数为负进而降低CMFD矩阵对角占优性，导致收敛不稳定或收敛发散的现象。基于广义等价理论的CMFD方法[6-7](generalized equivalence theory-CMFD, gCMFD)通过定义节块不连续因子与扩散系数修正因子，保证了CMFD矩阵对角占优性进而解决传统CMFD中可能出现的收敛不稳定问题。球床反应堆有极强的非均匀性，为避免在球床堆出现CMFD收敛困难或发散的现象，将gCMFD作为3D-MOC的加速算法,并使用广义最小残差算法[8-9]求解gCMFD方程，开发了gCMFD加速模块与MOCP程序耦合。

球床堆中存在大量随机分布的燃料球和冷却剂，对粗网格获取均匀化信息造成极大困难。为解决该问题，本文将燃料球细网格与粗网格进行距离判断以获得粗网内的平源区，同时针对MOCP中使用的CSG建模方法获得每个粗网内的冷却剂流平源区。最终得到gCMFD线性系统所需均匀化截面等参数。在实现粗网与平源区严格对应的过程中，可能出现平源区被分割的现象，文献[10]程序OpenMOC针对该问题有相应研究。

本文在数值验证部分，采用简单轻水堆基准题KUCA[11]对gCMFD的正确性与加速效果进行验证。自编直角几何球床堆芯和圆柱几何球床堆芯，初步验证gCMFD方法在不同几何下球床处理的正确性；高阶计算与低阶计算的等价性与加速效果。

1 gCMFD方程组与3D-MOC算法

gCMFD方法为3D-MOC提供更好的初值即中子通量密度和有效增殖因子以减少输运计算的迭代次数；同时，3D-MOC将平源区信息传递给gCMFD系统计算其均匀化截面。

1.1 三维特征线方法

将稳态下的中子输运方程改写为特征线形式，同时将梯度算符写为沿特征线的方向导数为：

(1)

源项为：

(2)

式中：E、Ω分别为能量、角度变量；χ、Σs、Σf分别为裂变能谱、散射截面与裂变截面。

将几何模型沿3个方向进行划分得到每个平源区网格，在能量、角度和材料方面引入不同的近似后，每个平源区网格的标通量表示为：

(3)

(4)

在MOC中，角通量为解析求解为：

Ψi,g,m,k(s1)=Ψi,g,m,k(s0)exp(-τi,g,m,k)+

(5)

式中：s1、s0为该平源区网格内特征线段的起点和终点；τi,g,m,k为光学厚度。

当3D-MOC求解完毕，通过平源区信息获得每个粗网格的均匀化截面并更新gCMFD线性系统特征值为：

(6)

1.2 圆柱几何下gCMFD

CMFD方法的本质是求解中子扩散方程。在圆柱几何下，对空间进行离散得到的网格如图1所示。

图1 圆柱几何空间网格Fig.1 Cylindrical meshes

gCMFD中，2节块间的中子流计算式为：

(7)

(8)

(9)

式中：f为节块不连续因子；g为扩散系数修正因子；h为相邻节块中心之间的距离；α为边界处的边界条件。真空边界条件时，α=0.5；全反射边界条件时，α=0。

如图1所示，圆柱几何下沿径向的粗网格会呈现逐渐变大的趋势，进而导致网格间的交界面积与距离发生变化，因此不同方向的交界面积与距离计算式分别为：

径向：

(10)

(11)

周向：

Ac,E=ΔrcΔzc,hc,E=r′cΔθc,hE,c=r′cΔθE

(12)

Ac,W=ΔrcΔzc,hc,W=r′cΔθc,hW,c=r′cΔθW

(13)

轴向：

Ac,B=r′cΔrcΔθc,hc,B=Δzc,hB,c=ΔzB

(14)

Ac,T=r′cΔrcΔθc,hc,T=Δzc,hT,c=ΔzT

(15)

其中，r′c表示该网格的几何中心：

(16)

几何中心表示某一规则网格最中心的位置，并通过几何中心的定义准确获得圆柱几何下离散网格间的长度，如图2。ri表示沿径向第i个网格边界中心的半径值。

图2 扇形网格几何中心示意Fig.2 Geometry center of shell mesh

同时，周向网格在边界处呈周期性边界条件，即：周向的起点网格与终点网格相邻，此时网格的耦合扩散系数用式(8)计算；径向第1层的网格，沿“SOUTH”方向没有交界面积，当前网格沿“SOUTH”方向的耦合扩散系数为0。

在gCMFD方法中强制高阶系统的中子流等于低阶扩散系统的中子流并以此计算NDF与MDF[6，12]。根据算法1在MOCP中加入新的功能以判断每条特征线段穿过粗网的方向。

算法1：find Boundry Segments CrossCM

Intersections←findIntersection(segment,CM)

ifIntersections is 1then

forallPlane∈CMdo

Length1←minDistancePointToSurface(segment,Plane)

endfor

elseifIntersections is 2then

forallQuadratic∈CMdo

Length2←minDistancePointToSurface(segment,Quadratic)

endfor

minLength←min(Length1,Length2)

Direction←setSegmentDirection(minLength,segment)

end

获得每个粗网格的节块不连续因子与扩散系数修正因子后可根据式(8)、(9)计算耦合扩散系数，最终得到gCMFD的线性系统：

(17)

2 球床模型处理

在球床堆中实现gCMFD加速的关键在于保证平源区网格与粗网的严格对应。球床堆没有类似轻水堆中栅元的结构，而在球床内存在大量随机分布的燃料球及冷却剂流，因此在球床堆获得粗网格内平源区信息是极为困难的。由于MOCP使用CSG建模方法的优势，各个几何所具有的面包含该几何相应的关键信息，为实现球床堆内的gCMFD加速提供了可能。

2.1 构造实体几何

构造实体几何是一种通过简单模型的“布尔”运算来构造复杂几何体的常见实体建模方法。在MOC的特征线追踪过程中计算点、线、面等所需要的参数信息均可通过CSG构造的复杂几何体所包含的面得到，因此CSG非常适合作为MOC的几何建模方法。

CSG方法建模的过程为：1)定义平面并表示相应半空间；2)对半空间进行“布尔”运算构造简单几何；3)将简单几何进行“布尔”运算构造复杂几何。

在三维空间中，某一平面或曲面将三维空间分为2部分：其中F+表示该面所定义出的正半空间，F-表示该面定义出的负半空间，如图3所示。

图3 三维面与半空间构成Fig.3 Construct of 3D surface and half space

通过对半空间进行“布尔”运算可构成不同的简单几何体；对简单几何进行“布尔”运算即可获得所需要的复杂几何。需要注意的是，为了保证建模的正确性，尽量通过不同的简单几何来构成复杂几何，并尽可能减少“布尔”运算的次数。在CSG方法中，复杂几何的最基本单位是各个不同的面，MOCP程序针对球床堆芯的复杂模型，提供不同类型的面，进而提高了描述几何的能力。CSG方法构造复杂几何的思路如图4所示。

图4 复杂几何构成Fig.4 Construct of complex geometry

通过CSG，球床堆芯内的燃料球由不同球面定义出的负半空间构成；堆芯容器由不同面方程定义出的半空间经过“布尔”运算构成；堆芯内的冷却剂流由堆芯容器与球面定义出的正半空间取“交”运算组成。半空间均保存相应模型的几何信息，如：球心位置、燃料球半径和平面参数，为粗网格获得平源区信息打下基础。

2.2 燃料球及球床间隙处理

在轻水堆中，粗网格主要是一个个栅元，对每个栅元进行划分后得到粗网格内的平源区。而在圆柱几何球床堆中，只能对整体圆柱几何按照轴向、周向和径向进行划分，得到虚拟的粗网格。因此，如何获得虚拟粗网格中的平源区及相应信息就是实现球床堆内gCMFD加速的关键。

在MOCP中，燃料球的处理方法通常是将燃料球沿径向划分后再进行卦限划分得到形状为“BALL”或“SHELL”的平源区网格。根据这些网格的面参数与粗网进行距离关系判断，若为“相交”或“属于”关系，则认为该细网格或部分细网格属于当前粗网格的平源区，距离关系l为：

(18)

式中：A、B、C、D表示平面参数；x0、y0、z0表示球心坐标。

获得燃料球与粗网各个面的距离后，通过半空间判断燃料球是否与粗网存在交集，二维示意如图5所示。

编号为1和6的燃料球与粗网各个面的距离均小于半径，因此认为属于该粗网；编号为3的燃料球球心在粗网内，且与粗网距离小于半径，因此认为该燃料球与粗网“相交”并将部分粗网面赋予该燃料球；编号为4的燃料球球心在粗网外，由式(18)可知，二者距离小于燃料球半径，将粗网的圆柱面赋予燃料球。

图5 粗网与平源区对应Fig.5 Map FSR and course mesh

当所有粗网与平源区均遍历完毕，对粗网格内形状为“BALL”或“SHELL”平源区的外层半空间取反，并加上粗网格的半空间，组成该粗网格内的冷却剂流平源区，见算法2与算法3。最终每个粗网格均包含若干燃料球平源区与一个冷却剂流平源区。

算法2：map Fuel Balland Course Mesh

forallcm∈Course Mesh_mapdo

locationType←calcLenthFSRandCM(fsr, cm)

iflocationType is intersectthen

new_fsr←copySectionCell(fsr)

addFSRintoCM(new_fsr)

elseiflocationType is belongthen

new_fsr←copyCell(fsr)

addFSRintoCM(new_fsr)

end

endfor

经过以上处理，保证粗网与球床区域内细网的严格对应，但产生了许多不规则网格。不规则网格的体积通过特征线积分得到：

(19)

算法3：create Inter Space FSR

forallcm∈CMdo

interspace←createNewCell()

forallfsr∈fsrStoryInCMdo

addSurface(interspace,fsr)

endfor

addSurface(interspace,cm)

addFSRintoCM(interspace_fsr)

endfor

值得注意的是，SHELL类型的平源区由2个球面组成，且半空间相反，因此在判断距离时要考虑SHELL的空心部分；如图5所示，与粗网“相交”的燃料球被赋予新的半空间，原始平源区被分割成多个新平源区，此时需要将旧的平源区消除以保证特征线追踪模块的正确性。

3 加速模块数值验证

在数值验证部分，首先选择一个简单轻水堆模型用于验证gCMFD加速模块在MOCP中的加速效果与正确性。由于实际球床堆几何模型复杂，本文仅针对球床堆中常见的直角几何与圆柱几何建立对应的gCMFD线性系统。因此本文自编简易直角与圆柱几何球床堆模型初步验证gCMFD加速方法在球床堆中的可行性与加速效果。

3.1 KUCA基准题模型

KUCA是一个小型轻水堆的两群输运问题，其1/4堆芯结构如图6所示，具体截面信息见表1、2。该基准题有2个状态：1)无控制棒，且控制棒区域为真空；2)控制棒全部插入堆芯。本文选择第2个状态进行验证。

图6 KUCA基准题1/4堆芯结构Fig.6 Configuration for CUKA benchmark

第2种状态的KUCA基准题全堆尺寸为50.0 cm×50.0 cm×50.0 cm，布置特征线时选择Lee求积组(4阶，共24个方向)。划分平源区时，每个细网格宽0.5 cm，3D-MOC的收敛标准为10-5。在加入gCMFD加速模块时，对该基准题的3个方向进行了多种划分形式，gCMFD的收敛标准为10-8。

表1 KUCA基准题截面信息Table 1 Cross section information of KUCA benchmark

表2 KUCA基准题散射截面信息Table 2 Scattering cross section info of KUCA benchmark

在验证过程中使用蒙特卡罗程序对相同模型进行计算作为KUCA基准题的参考解，特征值为1.038 93，耦合gCMFD加速模块的MOCP程序详细计算数据见表3、4。将gCMFD按不同划分方式进行计算，所得特征值均与高阶MOC计算保持一致，为1.038 99。可以看出，gCMFD实现与高阶计算严格等价的同时，极大降低3D-MOC计算的计算代价。

表3 KUCA基准题计算特征值数据Table 3 Eigenvalues of KUCA benchmark

表4 归一化通量分布Table 4 Distribution of normalized flux

3.2 自编球床堆模型

为验证gCMFD加速算法在不同几何下的准确性与加速效果，同时对算法2及算法3在球床部分的正确性进行验证，自定义直角与圆柱几何下的球床堆芯模型。

3.2.1 直角几何球床堆芯

该基准题为小型直角几何球床堆的两群输运问题，堆芯结构如图7所示，除中心的球床部分，还有外围的反射层，具体截面信息见表5、6。

图7 直角几何球床堆芯模型Fig.7 Cuboid pebble-bed reactor geometry

表5 直角几何球床堆模型截面信息

表6 直角几何球床堆模型散射截面信息

该基准题全堆尺寸为30.0 cm×30.0 cm×30.0 cm，球床尺寸为15 cm×15 cm×15 cm，球床内共172个燃料球，燃料球半径为2.99 cm，球床内的间隙为“void”，即冷却剂材料。在布置特征线时选用Lee求积组(4阶，共24个方向)。每个燃料球沿径向划分6层，再进行卦限划分，即单个燃料球划分出48个平源区。3D-MOC收敛标准为10-5，gCMFD收敛标准为10-8,特征值为1.185 76。详细计算结果见表7。

表7 自写直角几何球床堆特征值计算数据Table 7 Eigenvalues of cuboid pebble-bed geometry

3.2.2 圆柱几何球床堆芯

为验证圆柱几何下gCMFD系统的正确性，自编圆柱型球床堆芯，该基准题为两群输运问题。堆芯结构如图8所示，截面信息见表8。

圆柱堆芯半径为49.462 5 cm，高79.14 cm，球床部分圆柱半径为32.975 cm，高65.95 cm。球床内含有1 080个燃料球，燃料球半径为2.99 cm，燃料球间间隙为“void”。特征线布置中，每条特征线宽0.5 cm，燃料球划分方式同为沿经向划分6层后进行卦限划分。3D-MOC收敛表示为10-5，gCMFD收敛标准位10-8。特征值为1.260 56，其余计算结果见表8。

图8 圆柱几何自写球床堆模型Fig.8 Cylindrical pebble-bed reactor geometry

表8 圆柱几何球床堆特征值计算结果Table 8 Eigenvalues of cylindrical pebble-bed geometry

以上数值结果表明，gCMFD加速模块编写正确，且在保证不同分辨率系统严格等价的同时能够对3D-MOC输运计算起到良好的加速效果。

4 结论

1)本文中针对球床的处理方式正确，保证粗网与燃料球细网格严格对应，能够准确计算每个粗网均匀化截面等信息的同时，为细网格提供了更好的初值以加速收敛。

2)初步实现gCMFD加速方法在球床反应堆的应用，同时将gCMFD的几何限制拓展至圆柱几何，大幅提高了3D-MOC的计算效率。未来可针对圆柱几何特殊结构进一步改进gCMFD方法从而实现更快速准确求解。