去射线效应堆外探测器响应函数计算研究

徐龙飞，沈华韵，魏军侠，曹良志，郑友琦

(1.北京应用物理与计算数学研究所，北京 100094；2.西安交通大学，陕西 西安 710049)

反应堆堆外探测器可用于测量反应堆功率分布，探测器的响应函数表征了堆芯不同位置处单位强度裂变中子源对探测器读数的贡献。一般，堆外探测器响应函数可通过正向计算或共轭计算两种方式获得。正向计算的基本思想是：重复在堆芯不同位置处设置单位强度裂变中子源，并通过正向输运计算得到堆外探测器灵敏区的中子注量率，最后将不同位置处的计算结果进行归一化处理以得到探测器对各位置的响应函数[1]。正向计算方法的思想简单明了，但需对每个关心的区域单独进行输运计算，因此过程繁琐且计算量巨大。共轭计算方法则通过在探测器位置处设置共轭源项，进行1次共轭输运计算即可得到每个位置处的响应函数[2]。共轭算法由于计算量小，在工程实际中得到了广泛应用。中国核动力研究设计院[3]、上海核工程研究设计院[1]、清华大学[4]、上海交通大学[5]、西安交通大学[6]等单位均基于共轭输运算法对堆外探测器响应函数的计算展开了研究。目前，共轭输运计算主要采用SN方法和蒙特卡罗方法。SN方法具有计算速度快且可一次性获得全局共轭通量的优点，但射线效应[7]严重影响了SN方法的计算精度。西安交通大学郑友琦等[8]在利用二维综合SN方法计算堆外探测器响应函数时发现射线效应对不同组件响应函数带来的计算偏差可高达50%。蒙特卡罗方法具有良好的几何适应性且不存在射线效应问题，但堆外探测器响应函数计算属于深穿透问题，蒙特卡罗方法要获得可靠的计算结果需投入大量的粒子进行模拟，因此计算效率较低[9]。

基于以上分析可发现，克服SN方法在共轭计算中的射线效应问题对于改善堆外探测响应函数的计算精度和计算效率具有重要的工程实际意义。本文基于首次碰撞源射线效应消除方法[10]研究共轭输运计算中的射线效应消除策略，并针对首次碰撞源方法在工程实际中应用时所面临的计算量大、内存需求高等问题研究相应的并行化计算方法与动态内存管理方法。

1 堆外探测响应函数的共轭计算方法

1.1 响应函数定义

假定某堆外探测器灵敏区的中心位置为r0，灵敏区的有效体积为Vd，灵敏区对中子的响应截面为σd，那么探测器灵敏区的响应Rd可表示成：

(1)

现考虑在堆芯活性区ri位置处有一各向同性的单位强度裂变中子源，由该位置处裂变产生的中子到达堆外探测器灵敏区的中子通量密度记为ψi(r,E,Ω)，那么由ri位置处裂变产生的中子对探测器响应的贡献可写成：

Rd(ri→r0)=

(2)

一般，探测器的读数是探测器对堆芯功率分布的响应，即：

(3)

其中：P(r)为堆芯三维功率分布；V为堆芯活性区的体积；W(r)为探测器对功率的空间响应函数。根据反应堆功率与中子通量密度之间的关系，探测器的空间响应函数可表示为：

(4)

1.2 共轭输运方程

为简明表达，首先定义如下内积算符：

(5)

现以σd(r,E)为共轭源，则共轭输运方程可写成如下形式：

L*ψ*(r,E,Ω)=σd(r,E)δ(r-r0)

(6)

其中，L*为共轭输运算子。式(2)中的ψi(r,E,Ω)需满足如下正向输运方程：

(7)

其中：δ为狄拉克算符；χ(E)为裂变谱。根据共轭输运方程的性质[11]：

〈ψ*Lψ〉=〈ψL*ψ*〉

(8)

可推导出：

(9)

其中，φ*为共轭标通量。因此，只需知道堆芯的共轭标通量分布和裂变谱并结合式(4)与式(9)便可计算出探测器的响应函数。

2 共轭输运射线效应消除理论

2.1 共轭首次碰撞源方法

系统中任一位置处的粒子均可看成由两部分组成：一部分是由源项发射出来未经碰撞而直接到达该位置处的；另一部分则是由非源项位置散射到该位置处的[12]。依据这个原则，共轭通量可表示成：

(10)

将式(10)代入共轭输运方程中可得：

(11)

(12)

(13)

(14)

由于式(12)不包含散射源项，因此可解析求解。当共轭外源为各向同性时，式(12)的解析解可表示为：

(15)

(16)

其中：rp为未碰撞共轭源所在位置；-Ω为从位置r到rp的方向向量。将式(15)求得的未碰撞共轭通量代入式(13)，并作为已知共轭散射源项，即可求解出碰撞共轭通量，最后将两部分相加即可得到真实共轭通量。

2.2 快速径迹追踪方法

在求解未碰撞共轭通量时，需计算总截面沿追踪路径的积分τ。而τ的计算需对每条从点源到网格点的射线进行追踪，求出射线依次与每个网格的相交长度。为快速求解τ，本文研究了射线与网格相交的几何关系(图1)。对于方向为Ω(μ,η,ξ)的某射线，它必从网格的某个面穿出，穿出面取决于网格宽度在射线上投影长度的大小。当射线从X面穿出时，则必有|OA|<|OB|成立，即Δx/μ<δz/ξ(图1a)。同理，当射线从Y面穿出时，则必有Δy/η<Δz/ξ成立(图1b)。而当射线从Z面穿出时，有Δz/ξ<Δx/μ成立(图1c)。因此，只需求出(Δx/μ,Δy/η,Δz/ξ)中的最小值就可判断射线将从哪个网格面穿出，而射线与网格的相交弦长也可通过网格各方向的宽度在射线上的投影求得。

2.3 并行化计算方法与动态内存管理

式(15)与(16)的求解是相互独立的，因此具有天生的并行性。采用空间区域分解的方式对未碰撞共轭通量的求解进行并行化。以图2所示的二维网格为例，在左下角网格中心处放置1个共轭外源，现假定对该二维网格做2×2空间区域分解。那么，每个核只需计算各自子区内对应网格角点的未碰撞共轭角通量，同时每个核仍需保存全局网格的材料号以及所有材料的总截面。图2中的实网格即为每个核对应的子区域，虚网格即为每个核需额外存储的网格材料号和每种材料的总截面。本文采用MPI[13]对首次碰撞源进行并行计算，当式(15)求解结束时，存储全局网格材料号和每种材料总截面所占用的内存立即释放。

图1 射线与网格相交的几何关系Fig.1 Geometry relationship between tracer and computational mesh

图2 并行计算区域分解示意图Fig.2 Spatial domain decomposition in parallel calculation

由于τ与点源数、网格数以及能群数相关，其所需要的存储是巨大的。假设某一问题的网格数为6×106、点源数为50、能群数为50，τ的类型为双精度浮点数，那么其所需的内存约为：

6×106×50×50×8byte≈112G

(17)

若不进行射线效应消除，即便散射源的展开阶数为P3，存储通量矩所需的内存也仅为：

6×106×50×(3+1)2×8byte≈36G

(18)

因此，进行射线效应修正所需的计算内存至少为148G。随着外源数量的增多，所需内存将线性增加。个人计算机和中、小型工作站通常难以满足这样的内存需求。观察式(15)可发现，每个τ在计算完对应的未碰撞共轭通量后就不会再被使用。基于这一特性，本文在并行计算时计算出每个核存储τ所需的内存，当每个核的τ内存需求大于1G时，便以二进制格式将其并行写入磁盘中，待到需使用时再以并行方式读入到内存。由于每个τ只需读写1次，因此不会消耗大量的时间。基于这样的内存管理，程序将具备较好的鲁棒性。

3 Kori-1压水堆探测器响应函数计算

3.1 模型简介

Kori-1是由西屋公司设计的压水堆[14]，共含有121盒燃料组件，组件宽度约19.41 cm，堆芯活性区的直径和高度分别约250.00 cm和365.76 cm。Kori-1的整体结构[6]如图3所示，4对BF3探测器分别布置在1/8对称线上的堆腔处，每对探测器由上、下两个完全相同的圆柱形探测器组成，探测器中轴线距堆芯中心点的距离约211.93 cm，探测器半径约3.90 cm，单根探测器的高度约150.49 cm，每根探测器的外壁上均涂有0.07 cm厚的铝层。

图3 Kori-1的结构示意图Fig.3 Configuration of Kori-1

由于该堆沿x轴和y轴对称，因此只需计算1/4堆芯即可。基于文献[8]采用二维SN程序DORT[15]对该问题所作的角度空间网格敏感性分析的结果，本文在计算时采用了S16全对称求积组，空间网格划分为174×174×192。文献[8]将BUGLE-96数据库[16]的p1截面转化为MCNP的多群截面，并将整个堆芯活性区沿轴向分成16层，利用多群MCNP统计了每层上燃料组件的共轭通量密度，同时给出了第2、5、9层的堆外探测器响应函数。为使计算结果具有可比性，本文也采用了BUGLE-96数据库的p1截面，同时给出了相应燃料层的响应函数计算结果。

3.2 结果分析

本文首先以1对上、下对称的堆外探测器为共轭源，分别采用三维SN方法(SN)、首次碰撞源修正SN方法(SN-FC)以及多群蒙特卡罗方法(MCNP)对其进行共轭输运计算，并比较图4a中所示的A、B以及C组件的轴向共轭通量密度分布，结果分别如图5、6、7所示。当以1对探测器为共轭源时，堆芯中平面的共轭通量密度最大，此时SN、SN-FC的计算结果与MCNP的偏差也达到最大。SN-FC的计算结果相对于SN的计算结果均有所改善，距探测器越远的组件改善效果越明显。

图4 组件排布示意图(a)及单个探测器下由MCNP程序计算所得的各层燃料区域的归一化响应函数值(b)Fig.4 Lattices distribution (a) and normalized response function value obtained by MCNP code (b)

图5 A组件轴向共轭通量密度分布Fig.5 Axial adjoint flux density distribution of lattice A

图6 B组件轴向共轭通量密度分布Fig.6 Axial adjoint flux density distribution of lattice B

图7 C组件轴向共轭通量密度分布Fig.7 Axial adjoint flux density distribution of lattice C

以上部探测器为共轭源，SN-FC和SN计算得到的第2、5、9层燃料区(距堆芯活性区中平面的高度分别为148.59、80.01、11.43 cm)归一化响应函数值如图8所示。以文献[8]的计算结果(图4b)为参考解，图9示出了SN-FC和SN两种方法在第2、5、9层燃料区中归一化响应函数值大于0.000 01的各区计算偏差。对于图4a中红色外围组件，SN-FC和SN的最大计算偏差分别为5.21%与19.91%，SN-FC的计算精度有了明显的提高。共轭SN方法的误差主要来源于射线效应和输运方程的离散化，SN-FC方法能有效消除射线效应所带来的误差。同时，SN-FC方法在求解未碰撞通量时需将体源转化为点源，这一过程也会引入相应的误差。但堆芯活性区距堆外探测器较远，点源近似引入的误差几乎可忽略。另外，由于蒙特卡罗方法的统计涨落，在对称排布的组件上，MCNP的计算结果并不完全相等，而SN-FC和SN的计算结果在对称组件上完全相等。因此，图9中每一燃料区的计算偏差均是对称燃料区取代数平均后的结果。计算过程中，SN-FC和SN均采用72核并行，计算时间分别为20.65 min和14.07 min。MCNP在36核下的计算时间为2 731.02 min，SN-FC的计算效率约为MCNP的66倍。

图8 SN-FC与SN计算所得的归一化响应函数值Fig.8 Normalization response function value obtained by SN-FC and SN

图9 归一化响应函数值大于0.000 01的各燃料区的计算偏差Fig.9 Calculation deviation of fuel regions with normalization response function value of higher than 0.000 01

图10示出了A组件各群共轭通量的MCNP统计误差，从图10可看出，随着能量的降低，MCNP的统计误差越来越大。图11示出了A组件各群共轭通量占总共轭通量的比值，随着能量的降低，各群比值近似呈指数衰减。这主要是因为中子的能量越低，越难以到达堆外探测器灵敏区。另外，前10群的共轭通量之和占据总共轭通量的95.60%，而A组件前10群的共轭通量在第1层和第8层的统计误差均在3%以内，因此，可认为MCNP的计算结果是可靠的。SN-FC对A组件的响应函数的计算偏差略大于MCNP的统计误差，而SN方法对A组件响应函数的计算偏差则远大于MCNP的统计误差，进一步表明射线效应严重影响堆外探测器响应函数的计算结果。

图10 MCNP程序对A组件各群共轭通量的统计误差Fig.10 MCNP code’s statistical fluctuations of each group adjoint flux for lattice A

图11 A组件各群共轭通量占总共轭通量的比值Fig.11 Ratios of each group adjoint flux over total adjoint flux

4 结论

本文基于共轭首次碰撞源射线效应修正的三维SN方法计算了韩国Kori-1压水堆的堆外探测器响应函数，计算结果表明：射线效应严重降低了传统三维SN方法对堆芯活性区响应函数的计算精度。以MCNP的计算结果为参考解，SN-FC方法相比于传统SN方法有效地改善了压水堆堆外探测器响应函数的计算结果，SN-FC方法对外围组件探测器响应函数的计算偏差均在6%以内，同时SN-FC方法的计算效率相对于MCNP提高了66倍。SN-FC方法对于工程实际中的堆外探测器响应函数计算具有重要意义。此外，本文提出的动态内存管理方法与并行化算法，有效解决了首次碰撞源方法在大规模问题中的内存需求与计算效率问题，相比于未做射线效应修正的SN方法，计算效率只降低了46%。