液态金属冷却快堆堆芯物理分析软件LoongSARAX的验证与确认

路 瑶,杜夏楠,李爱鑫,高杰豪,陈文杰,郑友琦,王永平,吴宏春

(西安交通大学 能源与动力工程学院,陕西 西安 710049)

在第四代反应堆中,快堆凭借高的平均中子通量具有增殖裂变核燃料和嬗变长衰变期锕系元素两大特点和优势[1-2],是我国核能三步走战略的重要环节。与传统的核电厂热堆相比,快堆通常设计得较为紧凑,以提高其核燃料增殖性能,因堆芯内功率密度高,需要采用换热能力足够强的冷却剂带走堆芯的热量。目前在建或在运行的快堆均采用液态金属钠或铅作为冷却剂,这得益于液态金属强大的载热能力。同时采用液态金属冷却的快堆设计容易实现小型化,在海、陆、空、天以及特种同位素生产等领域均有迫切需求。

为提高堆芯性能,液态金属冷却快堆设计与传统压水堆相比,在燃料形式、能谱复杂程度以及堆芯尺寸规模等方面都具有非常显著的差别,其表现出的特征使得适用于传统压水堆设计开发的堆芯设计软件和方法不再适用。因此,快堆软件的自主化研发对于我国实现快堆稳步发展、技术弯道超车具有重要意义。

1 核电软件验证与确认方法

世界上很多国家有不同的核电软件,相对于国内,国外的核电起步较早,不同的大型核电技术公司有针对其核电软件独特且完善的验证确认方法[3],如美国西屋公司将核电软件的验证与确认称为合格性测试,法国阿海珐公司的核电软件——堆芯燃料管理软件分为两部分,分别是组件程序APOLLO2和堆芯程序ARTEMIS,分别用数值验证和基准题的方式进行验证。

在核工业领域,有6种常用的软件确认方法[4],分别是解析方法(BASICS)、程序对标(CODE)、电厂测量数据比对(PLANT)、实验测量(TEST)、审查测试(LICENSING)和回归测试(REGRESSION), 针对以上6种方法的使用频率做统计,可知审查测试、电厂测量数据比对与实验测量这3种确认方式应用相较广泛[5]。因此,本文使用实验测量的方式对LoongSARAX进行确认。

2 LoongSARAX的验证与确认策略

2.1 LoongSARAX程序简介

NECP-SARAX是西安交通大学核工程计算物理实验室自主开发的先进反应堆中子学分析计算系统[6-8],LoongSARAX程序则是在NECP-SARAX程序的基础上,针对液态金属冷却快堆的设计计算与安全审评需求,定制化开发的版本。该版本适用于装载二氧化铀、MOX以及金属燃料,采用液态金属钠、铅或铅铋作为冷却剂的快中子反应堆堆芯物理分析计算,主要包括截面产生模块TULIP[9]和堆芯稳态分析模块LAVENDER[10]。

2.2 LoongSARAX程序的验证与确认方法

LoongSARAX软件包括组件程序和堆芯计算程序,若将软件称为系统,则组件程序和堆芯计算程序称为子系统,组件程序包括组件共振计算模块、组件输运计算模块等,堆芯计算程序包括中子输运计算模块和反应性系数计算模块。

LoongSARAX的验证流程如下:将整个程序分成不同的模块,确定每个模块验证所需要的基准题并生成验证矩阵,在此基础上,再分为3个阶段分别进行模块验证、子系统验证和系统确认,按照验证计划完成验证。LoongSARAX子系统各模块及其验证方法列于表1。

表1 LoongSARAX子系统验证Table 1 LoongSARAX subsystem verification

本文对子系统中每个模块进行单独的验证后,将所有模块集成一个完整的计算程序,并对其进行软件的系统确认。本文采用公开文献以及评估反应堆物理基准实验的国际手册中的相关零功率实验对LoongSARAX程序进行系统确认,其中临界反应性、能谱特征、反应性效应、反应性系数、动力学参数以及反应率分布所采用的确认方法均为实验测量。

本文在确认子系统验证和系统确认的具体内容后,搜集每个模块所需基准题,并建立验证矩阵和确认矩阵,如表2、3所列。

表2 LoongSARAX验证矩阵Table 2 LoongSARAX verification matrix

表3 LoongSARAX 确认矩阵Table 3 LoongSARAX validation matrix

3 LoongSARAX的验证、确认及不确定度量化

3.1 组件程序验证

TULIP组件程序主要验证3个模块,分别为共振计算模块、输运计算模块、少群截面均匀计算模块[11]。3个模块的验证均是通过对例题进行校算来实现的,例题的参考解kinf、能群截面及通量由MCNP程序计算产生。为验证共振计算模块和少群截面均匀计算模块,利用MCNP程序统计1 968群以及均匀化33群的核素截面,并展开对比分析。

本文在验证时选取了不同液态金属冷却剂(钠、铅)、不同富集度情况下均匀模型组件、非均匀模型的一维平板组件及六边形组件。分别用TULIP程序和MCNP程序对均匀模型下各燃料组件进行计算并进行对比,包括kinf、总截面及裂变截面(235U、238U、239Pu、240Pu)、冷却剂Na或Pb的总截面、56Fe的总截面以及归一化能谱,截面及能谱都包括1 968群和33群。以ZPPR17A燃料组件均匀模型组件为例,部分截面计算结果如图1所示,能谱计算结果如图2所示。均匀模型所有组件的kinf计算结果如表4所列。C/E为计算值与实验值的比值。

图1 均匀燃料组件截面计算结果Fig.1 Calculation result of cross section for homogeneous fuel assembly

图2 均匀燃料组件能谱计算结果Fig.2 Calculation result of energy spectrum for homogeneous fuel assembly

表4 均匀燃料组件kinf计算结果Table 4 kinf calculation result of homogeneous fuel assembly

从图1可看出,中高能量段时两个程序的计算结果吻合较好、偏差较小,低能量段时偏差较大。从图2可发现,低能量段时通量接近于0,这是导致粒子输运程序统计偏差较大的原因,因此低能量段的结果不作为参考。由于MCNP和LoongSARAX在处理高能区阈能反应时的方式不同,因此导致图2在高能群位置相对偏差开始下降。其余核素的1 968群截面计算结果偏差的均方根在[0,2.1%]范围内,33群截面计算结果偏差的均方根在[0,7%]范围内。

在计算一维非均匀组件时,若是六边形组件,则采用等效一维圆环模型,即根据体积内核子密度守恒,将不同材料等效成半径不同的一维圆环在径向上进行堆叠,如图3所示,以从内到外第3圈为例,在面积不变的情况下,根据圆环面积的计算公式得出黄色区域的半径。以钠冷快堆MET1000一维模型ZPPR17A燃料组件为例,仅列出240Pu、23Na总截面计算结果,如图4所示。能谱计算结果如图5所示,一维模型所有组件的kinf列于表5。

图3 等效一维圆环模型Fig.3 Equivalent one-dimensional ring model

图4 非均匀组件截面计算结果Fig.4 Calculation results of cross section for heterogeneous fuel assembly

图5 非均匀组件能谱计算结果Fig.5 Calculation results of energy spectrum for heterogeneous fuel assembly

表5 非均匀燃料组件kinf计算结果Table 5 kinf calculation result of heterogeneous fuel assembly

由图4可看出,中高能量段TULIP程序与MCNP程序计算截面吻合程度较高、偏差较小,低能量段偏差过大。从图5可看出,由于典型快谱的低能量段的中子通量非常低,导致低能量段的统计偏差很大,因此,低能量段的结果不作为参考解。

3.2 堆芯程序验证

针对堆芯程序的验证,重点关注其中子输运模块的计算精度,为此,本文将基于TAKEDA基准题对堆芯程序中的中子输运模块展开验证。

以TAKEDA基准题2为例,采用程序对标的方法验证输运计算模块。TAKEDA基准题2模拟的是一小型快中子增殖反应堆(FBR)[12]。小型快中子增殖反应堆含有4个材料区域,分别为燃料区、轴向增殖区、径向增殖区和控制棒区。当控制棒提出时,采用钠来填充空出的区域。

计算时考虑两种情况:控制棒全提(情况1)和控制棒插入一半(情况2)。有效增殖因数及控制棒价值计算结果列于表6,平均中子通量密度偏差计算结果示于图6,其中G表示能群。可看出,LAVENDER与MCNP程序所提供的参考值吻合很好,keff偏差小于40 pcm,控制棒价值偏差小于2%,各区域、能群的中子通量密度偏差均小于5%。

图6 平均中子通量密度偏差计算结果Fig.6 Calculation bias results of average neutron flux density

表6 有效增殖因数及控制棒价值计算结果Table 6 Calculation results of keff and control rod value

3.3 LoongSARAX程序验证

以OECD循环快堆和超凤凰堆基准题为主要计算对象,在输运、功率、燃耗、反应性效应等计算功能方面进行LoongSARAX程序验证。

3.3.1OECD循环快堆 钚循环物理工作组(WPPR)指定了两种快速燃烧器基准设计(1种氧化物和1种金属)。本文通过计算金属燃料反应堆进行验证。金属燃料反应堆堆芯的功率是1 575 MWth,周期长度为365 d,容量系数为85%。反应堆包含420个燃料组件和30个控制组件,并排列成一堆芯为1/6对称的配置,燃料组件由铀、钚组成[13]。

本文分为5个燃耗步且每个燃耗步时长为73 d,进行1个周期的循环。计算结果列于表7,其中BOL表示寿期初燃耗重核质量,EOL表示寿期末燃耗重核质量。从表7可知,LoongSARAX的计算结果和其他单位的计算结果相比,偏差较小。

表7 燃耗重核质量及其变化计算结果Table 7 Calculation results of burnup heavy core mass and its variation

3.3.2超凤凰堆 超凤凰堆是世界上运行的第一座商业规模的快中子增殖堆,也是至今已建成的最大的钠冷却的快中子反应堆,该反应堆装载358个燃料组件,反应堆额定热功率为3 000 MW,堆芯燃料区分区装载有两种不同的MOX燃料,增殖区装载低富集度的二氧化铀增殖燃料[14]。LoongSARAX通过超凤凰堆计算验证了输运、功率及反应性效应。其中,参考解来自于基准题报告中Serpent的计算结果。

本文完成了13种堆芯状态下的反应堆堆芯临界计算。这13种堆芯状态分别在几何温度、截面温度、控制棒插入深度上存在差别。13种堆芯状态的稳态计算结果与参考解之间的偏差列于表8。

表8 堆芯稳态计算结果Table 8 Core steady state calculation results

从表8可看出,LoongSARAX的计算结果和参考解相比不存在较大的偏差。除堆芯状态较为复杂的算例9外,其他所有算例和参考解的偏差都在400 pcm以内。表明LoongSARAX程序在堆芯临界计算上存在较高的计算精度。

堆芯功率分布计算基于的堆芯状态为几何膨胀温度673 K,调节棒吸收体插入活性区高度的40%。计算结果详见文献[15],计算结果和参考解的相对偏差如图7所示。

图7 功率分布计算偏差Fig.7 Calculation bias of power distribution

如图7所示,功率分布计算的计算结果和参考解相比在活性区不存在明显的偏差。增殖区相对偏差较大,绝对偏差不大,可认为是由于增殖区功率较低导致的。因此,可证明LoongSARAX程序对功率分布的计算结果具有较高的计算精度。

本文将基于与堆芯临界状态相同的堆芯状态计算堆芯的反应性常数和反应性系数,其中包括燃料多普勒常数、钠密度系数、钠空泡系数、燃料轴向膨胀系数、包壳膨胀系数、组件壁膨胀系数、栅格膨胀系数。具体计算假设和计算方法参考基准题报告[14],反应性系数和反应性常数的计算结果和计算偏差如表9所列。

表9 反应性系数和反应性常数计算结果Table 9 Calculation results of reactivity coefficient and reactivity constant

整体来看,LoongSARAX的计算结果和参考结果不存在较大的偏差。因此,可证明LoongSARAX程序在计算反应性系数和反应性常数时存在较高的精度。

3.4 LoongSARAX程序确认

本文以JOYO以及ZPPR17A两个零功率实验装置为主要计算对象,对程序在确认过程中的对比参数进行介绍。

3.4.1JOYO反应堆 JOYO是日本的第一个实验反应快堆,并在1977年达到了初始临界[16]。JOYO的主要作用是改进快堆技术、进行辐照燃料,并创新将快堆投入实际使用的技术。JOYO是钠冷快堆,采用六边形组件,以铀钚混合氧化物为燃料。在稳定的温度和恒定钠流量条件下,该反应堆主要通过调整控制棒的位置来实现临界。当堆芯装载70个燃料组件时,简称为JOYO-70,其中有两根调节棒组件半插入堆芯。本文主要进行了临界反应性、控制棒价值、钠空泡反应性以及燃料替换反应性计算,结果列于表10～13。

表10 JOYO-70堆芯临界计算结果Table 10 JOYO-70 core critical calculation results

从表10可看出,SARAX对临界测量实验的建模计算与实验结果吻合较好,有效增殖因数偏差为365 pcm。

在本次计算中,采用反应性差值法进行计算控制棒价值。选取JOYO-70所有控制棒均为半插入状态时的βeff作为本次计算采用的βeff,所有控制棒都撤出(ARO)时的有效增殖因数的测量结果为1.024 95,再分别计算每根控制棒插入时的keff,与ARO状态的有效增殖因数作差,即可得控制棒价值,最终测量结果如表11所列。可看出,LoongSARAX计算控制棒价值的结果较好,偏差较小。

表11 JOYO-70控制棒价值计算结果Table 11 JOYO-70 control rod worth calculation results

通过将燃料组件替换为钠空泡实验组件的方式得到钠空泡反应性,如表12所列。可看出,钠空泡反应性数值的绝对值减小。由于[000]处于堆芯,燃料数量多,中子通量最多,由堆芯径向往外中子通量减少,引入负反应性减少,[4f1]、[6F1]均在燃料区外,处于增殖区,移除钠引起正反应性。实验值与测量值偏差较小,所计算偏差均小于实验不确定度。

表12 JOYO-70钠空泡反应性计算结果Table 12 JOYO-70 calculation results of sodium void reactivity

燃料替换反应性的所有测量都是在低功率运行(温度分布均匀)下进行的,平均冷却剂温度范围为240～250 ℃。计算结果列于表13,可看出,两者结果较为吻合。

表13 JOYO-70燃料替换反应性计算结果Table 13 JOYO-70 fuel replacement reactivity calculation results

3.4.2ZPPR17A反应堆 ZPPR17A是零功率物理反应堆[17],其特征是轴向不均匀,临界堆芯状态没有控制棒位置。它被认为是JUPITER-Ⅲ系列实验的参考堆型之一。ZPPR-17A的燃料是Pu-U-Mo板,钠为冷却剂,反射层材料为不锈钢。本文采用LoongSARAX计算的临界反应性为1.003 43,与实验参考值的偏差为333 pcm,本文还计算了控制棒价值以及反应率分布,结果如表14及图8所示。可看出,计算结果与参考值较为吻合。

a——239Pu在z=5.16 cm处沿x轴总反应率分布;b——238U在z=28.02 cm处沿x轴裂变反应率分布;c——235U在z=28 cm处沿y轴裂变反应率分布;d——238U在148-70位置处沿z轴俘获反应率分布图8 部分反应率分布计算结果Fig.8 Calculation results of partial reaction rate distribution

表14 堆芯控制棒价值计算结果Table 14 Core control rod worth calculation results

基准题报告中,在燃料组件或增殖组件中装入金属探测箔片,通过测量箔片的感生放射性活度,并结合已知的探测箔片位置和材料特性,可得到相应区域的反应率,包括轴向高度z=5.16 cm、z=28.02 cm处沿x轴反应率分布、计算组件位置位于148-50(堆芯中心位置)和148-70(外燃料区内,具体位置见文献[17])的z轴反应率分布以及轴向高度z=28 cm处沿y轴的反应率分布(包括总反应率、裂变反应率、俘获反应率)。利用LoongSARAX计算反应率,所得结果已归一化。图8中仅展示其部分反应率分布。可看出,反应率分布的计算值与实验测量值的偏差较小。

3.5 程序确认结果分析

在3.3节的基础上,LoongSARAX针对确认矩阵中的算例及其他钠冷快堆如中国实验快堆(CEFR)均开展了计算,结果表明,与实验测量值相比,堆芯有效增殖因数计算偏差均小于400 pcm,其他反应性与参考解的偏差均小于15%。受限于篇幅,其他计算结果可参考文献[18-19]。

特别地,针对CEFR的启动物理试验,在临界反应性、控制棒价值、钠空泡系数、等温温度系数以及燃料组件替代反应性测量过程中共产生了60余组临界状态。具体结果如图9所示。

图9 中国实验快堆各临界状态计算结果偏差Fig.9 Calculation bias results of CEFR core critical reactivity

基于以上结果,在不依赖于总体分布类型的情况下,本文应用非参数统计法对CEFR的临界反应性进行了不确定度量化[20]。

为了提高结果的准确性和可靠性,本文对偏差数据进行筛选,剔除可能由于失误或不可控因素导致偏差数据存在异常值或数据误差较大的情况。

如式(1)[21]所示,假设从具有连续累积分布函数F(x)的总体中抽取大小为n的样本,将样本值按增序排列为x1、x2、…、xn,对于包含在样本中从第r个最小值xr到第s个最大值xn-s+1之间的比例可表示为F(xn-s+1)-F(xr)。这个比例称为区间(xr,xn-s+1)的覆盖率。如式(2)[21]所示,计算这个覆盖率至少达到某个给定数值α的概率,α表示100%的总体被包含在区间xr和xn-s+1之间的概率,这个概率仅依赖于n和m。因此,在样本量n一定的情况下,筛选的具体个数m根据置信度、概率查表[21]获得。

un-m(1-u)m-1du

(1)

(2)

在99%的置信水平下,通过查表[21]可知,经剔除筛选偏差数据后,LoongSARAX对于钠冷快堆CEFR的计算偏差有90%的概率位于[-389 pcm,300 pcm]以内。

4 总结

针对液态金属冷却快堆的研发需求,本文基于LoongSARAX软件提出了堆芯核设计程序的验证与确认策略,通过搜集国际上关于液态金属冷却快堆物理计算基准题建立了相应的验证与确认矩阵,在针对LoongSARAX制定验证确认流程后,开展了组件程序、堆芯程序各个模块的验证以及LoongSARAX程序确认的计算。

本文主要获得了各堆芯及其在实验过程中不同状态点下的有效增殖因数,以及其他反应性的结果与偏差,并针对有效增殖因数开展了计算结果的不确定度量化。由计算分析可得,与实验测量值相比,LoongSARAX对于临界反应性、控制棒价值、钠空泡价值和反应率分布等结果的计算精度较好,可说明LoongSARAX在液态金属冷却快堆的计算分析中能够满足多种中子学参数的分析需求且具备较高的计算精度。