小型铅冷快堆堆芯物理计算软件的开发与临界实验验证

陈仁宗,周 琦,朱庆福,夏兆东,宁 通,马骁笛,孙 旭

(1.清华大学 能源环境经济研究所,北京 100084;2.中国原子能科学研究院,北京 102413)

铅冷快堆(LFR)是第4代核能系统论坛(GIF)选定的6种优选堆型之一[1],随着材料和设备技术的进步,LFR在中子物理、传热、安全等方面具有的独特优势不断凸显[2],已成为世界各国新型核能系统研发的重点[3]。当前LFR堆芯物理计算主要以确保精度的蒙特卡罗方法为主,进行全堆芯的建模计算,随着型号设计的不断深入,效率较低的蒙特卡罗全堆计算设计方法已无法满足LFR设计需要巨大计算需求。

如果参考压水堆堆芯物理计算的两步法[4],首先采用较为精确的中子输运理论给出栅元和组件的等效均匀化参数,再将其应用于基于中子扩散理论的堆芯计算,相较于目前常用的全堆芯蒙特卡罗精确建模计算,能够大幅提高计算效率。与压水堆或钠冷快堆(SFR)相比,小型LFR冷却剂中的铅、铋等核素对中子的弹性、非弹性散射更强,输运截面较大,但慢化能力更小,燃料元件栅距变化范围大,空间异性强,燃料组件的等效均匀化遇到的问题更加特殊[5]。

本文参考国内外研究成果,针对小型LFR特有的材料与能谱环境,利用开源的蒙特卡罗计算软件OpenMC计算多群截面,能够更好地描述栅元与组件在空间与能量上的局部特性。调用开源的有限体积法(FVM)软件OpenFOAM的各种求解器,开发中子扩散求解器DESOF,通过Python实现蒙特卡罗均匀化与扩散计算的耦合,形成完整的堆芯物理计算软件MCDESOF,并开展临界实验,利用实验数据验证MCDESOF的适用性。

1 蒙特卡罗方法生成多群截面

1.1 多群截面的蒙特卡罗统计方法

蒙特卡罗均匀化的核心内容是将求解输运问题得到的各种事件统计结果转化为均匀化常数,这些事件主要包括散射、裂变、俘获吸收、(n,xn)等反应,均匀化常数包括各种反应的群截面、输运修正、散射矩阵、各阶勒让德项,次级中子产生矩阵、裂变产额、裂变能谱,以及动力学计算需要的中子速度倒数、缓发中子先驱核产额与衰变常量等。为了简化表达式,引入内积算符〈·,·〉来代表能量、空间或角度的积分形式,这样第x类核反应的反应率计数估计可表示为:

(1)

式中:Σx为第x类核反应的宏观截面;φ为中子通量密度;r为中子的空间位置;E为中子能量;Ω为中子的运动方向;V为r的积分空间范围;S为Ω的积分空间角范围。

(2)

使用径迹长度法,通过设置反应率和中子通量密度的计数箱可计算得到式(2)的分子与分母,分别是k空间(栅元)、第x类核反应、第g能群(能量区间)的核反应率R和k空间、第g能群的中子通量密度φ。

1) 输运修正

式(2)考虑的散射是各向同性的,但实际散射存在各向异性,因此需要引入输运修正处理各向异性散射。首先,定义以勒让德多项式Pl(μ)展开的散射核为:

(3)

式中:Σsl为散射核;E′为散射前的能量;μ为散射角余弦。

对于k空间、Vk体积、g群能量范围勒让德多项式可定义为:

φ(r,E′,Ω)dE′dEdΩdr

(4)

式中,g′为散射前的能群。

(5)

式中,G为能群总数。

(6)

由于输运修正项必须通过模拟估计法计算得到,因此总碰撞反应率与中子通量密度也需要通过模拟估计法得到,确保计算结果准确性。

2) 散射矩阵

(7)

对于群间散射,不仅要考虑群间散射的概率,还要考虑中子散射前后的能量和角度分布,蒙特卡罗方法可通过设置能量与角度计数箱,统计群间散射率与散射后的角分布。参考Redmond[6]提出的计数箱设置方法,首先模拟用于统计群间散射的散射事件,根据入射粒子的能量通过散射事件计算出射粒子的能量、散射角余弦,其次入射及出射能量决定了哪个群间散射率计数箱储存信息,散射角余弦则决定了哪个群间散射角余弦接受信息。

类似输运修正(式(6)),散射矩阵的修正可表示为:

(8)

3) 需要统计的物理量

为了提供后续扩散计算所需要的扩散系数、移出截面、裂变截面等完整的输入参数,参考Boyd等[7]提出的计算思路实现。

1.2 利用OpenMC生成多群截面

利用OpenMC生成多群截面计算步骤如下。

1) 问题的预处理,主要目标是建立蒙特卡罗计算模型,根据问题所包含的几何、材料与均匀化计算的要求,根据所选蒙特卡罗软件的要求,利用Python按照特定的语法建立计算模型,准确定义几何、材料、截面、计算参数与计数箱的设置。

2) 蒙特卡罗计算,主要目标是生成多群截面所需的原始数据。调用蒙特卡罗软件进行计算后,将生成可直接读取的ASCII格式与便于处理的二进制格式的包含计算结果的文件,继续利用Python将其数据提取、处理并计算得到吸收、裂变、散射截面与散射矩阵、不同阶的勒让德散射项、裂变能谱与缓发中子参数等基础数据。

3) 多群截面的应用,主要目标是对基础多群数据进行处理,具体过程将根据扩散计算采用的近似处理方法进行处理,如进行输运修正,生成扩散系数、散射移出截面、散射移入截面等。具体的文件格式根据扩散求解器的构造利用Python进行适应性的调整。

1.3 计算流程

以OpenMC为例,本文利用蒙特卡罗方法计算多群截面并为多群中子扩散计算提供均匀化参数的计算流程,如图1所示。

图1 蒙特卡罗方法为中子扩散计算提供均匀化参数的流程

2 多群中子扩散软件开发

2.1 有限体积法与OpenFOAM

有限体积法又称为控制体积法,其基本思路与有限元方法类似,也是节点划分、方程离散、方程求解散步,不同在于将计算区域划分为一系列不重复的控制体积,并使每个网格点周围有1个控制体积,将待解的微分方程对每一个控制体积积分,得出1组离散方程,其物理意义就是因变量在有限大小的控制体积中的守恒原理[8]。

OpenFOAM[9]是一个完全面向对象的CFD类库,包含许多可执行文件,也称为应用程序包,从文件组织结构来说这些应用程序大体可以分为两大类:求解器和辅助工具。求解器用来求解连续介质力学中的某个特定问题,而辅助工具主要用来进行数据操作、辅助求解器完成计算任务。

2.2 DESOF开发

中子扩散方程可看作一个基础的偏微分方程,OpenFOAM针对中子扩散方程的求解通过式(9)的结构进行,主要包含拉普拉斯项和源项两类离散项,左边第1项泄漏项与第2项吸收项(实际上包含吸收与群件散射移出两部分)采用自带的隐性finiteVolumeMethod离散方式实现右边的源项neutronGroups包含群间散射移入项与裂变产生项。

-fvm∷Laplacian(D,φ)+fvm∷Sp(ΣA,φ)=

neutronGroups.Sp(neutron)

(9)

基于C++语言,根据OpenFOAM求解稳态中子扩散方程采用的代码形式,本文开发了基于OpenFOAM的扩散方程求解器DESOF。DESOF的基本计算流程如图2所示,主要分为3部分:OpenFOAM输入、网格生成及扩散方程的求解。

2.3 数值验证

选取二维IAEA基准题对DESOF进行验证。基准题是一简化的两群PWR基准题[10],它由241个组件、4种不同的材料构成。堆芯按双区布料方案布置,组件几何尺寸为20 cm×20 cm,堆芯外边界条件为真空边界。基准题1/4堆芯几何布置、DESOF构建的几何布置如图3所示。图3中,Jin为入射界面流。

图3 基准题(a)与DESOF(b)建立的几何模型

基准题keff的参考值为1.029 585,DESOF的计算值为1.029 400,二者偏差为18.5 pcm。将组件归一化功率分布计算值与参考值进行比较,如图4所示。由图4可见,最大相对偏差出现在边缘处的组件,为1.00%。从对比结果看,只要输入准确的均匀化参数,DESOF的多群扩散计算准确性能够得到保证。

图4 基准题的组件归一化功率分布

3 组件均匀化与验证

3.1 MCDESOF的集成与功能性检验

为了提高计算效率,本文基于Python开发了链接蒙特卡罗均匀化计算与多群扩散方程求解器的Pylink程序,实现便捷的蒙特卡罗均匀化建模、计算与数据处理,传递并进行扩散求解,形成了完整的MCDESOF软件。其中,蒙特卡罗均匀化建模模块将自动生成OpenMC运行所需的几何、材料、设置、图形和计数等部分的XML格式的输入文件,调用OpenMC完成蒙特卡罗输运计算;多群截面处理模块处理蒙特卡罗输运计算结果,将其处理运算得到多群扩散需要的均匀化参数;多群扩散计算模块对指定网格的几何排列方式逐一进行赋值、设定网格的边界条件以及图形化显示,并调用DESOF完成扩散求解;结果处理模块将扩算计算结果转化为keff和堆芯功率分布或组件分布,便于后续处理。

3.2 小型LFR燃料组件模型

参考各种LFR燃料组件的特点,本文构建了4种类型的组件,如图5所示。选取金属铀芯块进行计算验证,燃料富集度为15%,燃料包壳为不锈钢,冷却剂材质选择铅铋合金。燃料棒芯块外径为0.9 cm,包壳内径为0.92 cm,包壳外径为1 cm,燃料棒对边距为1.40 cm,燃料棒慢化剂外边界条件为全反射边界。

使用蒙特卡罗方法对4个计算案例进行输运计算,给出k∞与各群相对中子通量密度计算结果,作为本文提出的等效均匀化方法的对比参考值。蒙特卡罗计算考虑keff统计误差1σ小于10 pcm。采用3群能群结构(分界能分别为5.53×10-3MeV和0.8 MeV)与网距为1/2栅距的网格划分方式,利用MCDESOF进行了计算。蒙特卡罗均匀化计算考虑的k∞统计误差小于10 pcm,扩散计算以特征值偏差|ε|<1 pcm作为收敛判据。将MCDESOF的计算值与蒙特卡罗参考值对比分析,结果列于表1、2。

表1 MCDESOF计算值与参考值的k∞对比

表2 案例1 MCDESOF计算值与参考值的3群中子通量密度对比

由表1、2可见,k∞最大偏差超过400 pcm,最小偏差也超过了200 pcm,各群的归一化中子通量密度最大相对偏差超过1%,符合得不够理想,因此有必要对蒙特卡罗多群计算得到的均匀化参数的等效性进行核实,并进行适当修正。

3.3 超级均匀化方法的应用

为保证精度,均匀化过程必须保证一些重要的、反映堆芯物理特性的参量在均匀化过程中守恒,这些参量一般包括反应率、界面流和特征值。在实际情况中,全部满足3个参量守恒是不现实的,一般通过放宽约束条件或增加均匀化参数的自由度来达到或大致满足等效均匀化的要求,当前应用广泛的是广义均匀化理论(GET)和超级均匀化(SPH)方法。广义均匀化理论通过组件边界处的不连续因子保证界面流守恒,由Henry等[11]提出。Gunow等[12]提出不连续因子的具体实现方法,该方法目前在轻水堆工程计算中应用广泛。Koebke[13]提出超级均匀化方法直接调整截面保持反应率守恒。Hébert等[14]先后发展了三代超级均匀化方法。

相比广义均匀化理论,超级均匀化方法的优势在于无需额外的均匀化参数,这使得其在精细栅元计算中成为首选。本文使用超级均匀化方法,将中子扩散方程的右边简化为固定源项,在空间k、能群g中引入SPH因子μk,g,扩散方程通过修正移出项实现平衡,如式(10)所示。

(10)

式中:φg为g群中子通量密度;Σt,k,g为空间k处的g群总截面;Qk,g为空间k处的g群源项。

本文利用迭代算法得到SPH因子,首先将式(10)转化为迭代的形式:

(11)

(12)

针对案例1,采用超级均匀化方法对各群截面进行修正,在式(16)中设定ε1为50 pcm,ε2为0.5%。

用超级均匀化方法对4个案例的各能群均匀化参数进行修正,得到修正后各组件的k∞计算结果,并与蒙特卡罗参考值进行对比,结果列于表3。由表3可见,经过修正后k∞的偏差基本下降了1个数量级,最大偏差不超过50 pcm,充分显现出该方法对于确保均匀化等效性的作用。

表3 修正后MCDESOF计算值与参考值的k∞对比

4 临界实验验证

4.1 临界装置与实验方案

为了验证MCDESOF对于小型LFR堆芯物理计算的适用性,在中国原子能科学研究院铀棒栅临界实验装置的铅堆上开展了临界实验。该装置是为开展先进反应堆堆芯物理研究而构建的铅铋快堆零功率装置,陆续开展了铅基ADS次临界反应堆物理特性实验研究(启明星Ⅱ号)[15]和铅铋快堆基准性零功率实验研究[16],已成为获取小型LFR核心堆芯物理实验参数的重要实验平台[17]。

本次临界实验方案的堆芯中,以金属铀为核燃料,主要以金属铍为内部反射层,以石墨为外部反射层,以铅铋合金作为堆芯介质体,燃料棒、铍棒与控制棒等布置于打孔的铅铋介质体中。燃料元件按六角形栅格排布,4圈燃料构成1个组件,单个组件包含37根燃料棒,整个堆芯包含19个燃料组件。安全棒采用富集度为92%的碳化硼作为吸收体材料,布置在外围反射层组件中。临界实验堆芯布置如图6所示。

图6 铅堆临界实验堆芯布置

4.2 MCDESOF建模

1) 组件网格划分

在径向上,采用菱形网格划分方式,针对每个不同的组件,通过预先制定好的蒙特卡罗均匀化截面,将其填充到六角形的组件中。对于六角形组件的网格划分,由于基于OpenFOAM的MCDESOF本身具有完整的六面体网格划分方式,燃料组件和铍反射层组件的划分方式如图7所示。

图7 燃料组件(a)与粗铍棒组件(b)的网格划分

在轴向上,活性段划分100层,上下反射层各划分了20层水平方向的平行网格。

2) 蒙特卡罗均匀化与修正方案

根据堆芯中的各种组件的排列方式以及堆芯的对称性,组件可分为6种,如图8所示。其中:1号和2号组件,它们拥有较多数量的燃料元件,采用反射边界的组件计算模型,考虑SPH因子,为扩散计算提供修正后的多群截面;3号、4号组件,它们不含燃料元件,计算时将2号组件的1/2与3号、4号组件以及外围的石墨反射层组成“超组件”计算模型,内部采用反射边界,外部考虑泄漏边界,考虑SPH因子,为扩散计算提供修正后的多群截面;5、6号组件,控制棒不插入时,这两种组件等价于4号组件,不需要额外计算,一旦控制棒插入,采用的相同的方法,引入2号组件的1/2组成“超组件”计算模型,考虑SPH因子,为扩散计算提供修正后的多群截面。

图8 堆芯蒙特卡罗均匀化的区域设置

3) 堆芯扩散计算网格划分

对于堆芯扩散计算,考虑到对称性,采用菱形网格读取蒙特卡罗均匀化计算得到的堆芯各组件多群截面,进行全堆扩散计算。对于堆芯外围与石墨反射层接触的不规则结构的网格,如图9所示,以任意一处边界为例,首先计算出组件与边界处的交点1、3、5的坐标,然后使用三棱柱prism关键字对(1 2 3)以及(3 4 5)两个面所在的棱柱进行定义。由于1-3以及3-5的边分别为圆弧形,因此在进行网格定义的输入时,除点、面的定义外,还需增加对边的定义。

图9 堆芯边界处不规则网格的处理方式

4.3 临界实验结果与计算值对比

在启明星Ⅳ号上开展了临界装载实验和安全棒价值测量实验。临界装载实验堆芯共装载198根燃料元件,1个空孔道,其余位置被金属铍反射元件填充,此时反应堆倍增周期为82.9 s,反应性为47.1 pcm。安全棒价值测量实验采用了基于逆动态方法的落棒法[18]。落棒后使用逆动态方法得到1号安全棒价值为1 383.8 pcm。

利用MCDESOF对临界附近的反应性与安全棒价值进行了计算,同时利用OpenMC采用蒙特卡罗方法直接对全堆进行建模计算,核截面数据库选择ENDF/B-Ⅶ.0。计算的硬件平台为Dell T7920服务器,64核心并行计算,蒙特卡罗均匀化统计误差1σ为10 pcm。实验测量值、蒙特卡罗全堆建模计算值、MCDESOF计算值三者的对比列于表4。

表4 实验测量值、蒙特卡罗全堆建模计算值与MCDESOF计算值的对比

由表4可见,MCDESOF计算值、蒙特卡罗全堆建模计算值与实验测量值间的偏差相当,最大偏差不超过200 pcm。但是两种理论方法在安全棒、调节棒插入前的计算值大于实验值,而在安全棒或调节棒插入后计算值小于实验值,这是因为安全棒及其配套的结构比较复杂,需要提高建模的精度,特别是吸收体密度、填充范围与尺寸需要进一步核实,减少建模误差。

从计算时间上看,OpenMC采用蒙特卡罗方法直接对全堆进行建模计算需要约2 h完成,MCDESOF内部实施的组件蒙特卡罗均匀化计算,截面处理与扩散计算总共不超过25 min,计算时间小于蒙特卡罗方法直接计算的25%,计算效率有显著的提升。

5 结论

本文开展了小型LFR堆芯物理分析软件的开发与临界实验验证。利用开源蒙特卡罗软件OpenMC,实现栅元与组件的蒙特卡罗均匀化。基于开源有限体积法软件OpenFOAM开发了中子扩散求解器DESOF,通过Python实现蒙特卡罗均匀化与扩散计算的耦合,形成完整的堆芯物理计算软件MCDESOF,利用超级均匀化方法实现了组件的等效均匀化,通过4种LFR典型的燃料组件模型对MCDESOF进行了数值验证。在铀棒栅临界实验装置铅堆堆芯上开展了临界装载实验与控制棒价值测量实验,利用MCDESOF进行了针对性的建模,进行了全堆扩散计算,计算结果与测量结果进行对比,临界附近的反应性偏差小于100 pcm,安全棒价值相对偏差小于200 pcm,计算准确度达到蒙特卡罗全堆建模计算的水平,所需的计算时间则远小于蒙特卡罗全堆建模所需时间。

下一步将完善控制棒均匀化计算方法,引入误差传递公式,定量分析蒙特卡罗均匀化截面的不确定度对扩散计算的影响,以及开展三维时空相关的瞬态计算功能开发,支持小型LFR瞬态安全分析。