吴小飞,吴海成,刘 萍,葛智刚
(中国原子能科学研究院 核数据重点实验室,北京 102413)
燃耗计算中的核素积存量等响应参数是核反应堆设计中非常重要的积分参数,其准确度直接影响核装置设计的精度。燃耗计算响应参数的不确定度主要来源包括输入参数、物理模型和数值求解方法。随着核装置物理计算模型和计算机技术的发展,物理模型简化近似和数值求解方法引入的不确定度显著降低,输入参数引入的不确定度已经成为响应参数不确定度的主要来源[1]。核数据是燃耗计算最基本的输入参数,量化核数据不确定性引起的燃耗计算响应参数的不确定度,对提高反应堆设计和乏燃料后处理/贮存/运输的安全性、经济性有着重要作用。
目前国内外开展的燃耗计算不确定度分析方法主要分为微扰方法和随机抽样方法。基于燃耗微扰理论[2]的程序较少,主要有日本北海道大学的CBZ程序[3]、日本JAEA的MARBLE程序[4]等。基于随机抽样方法的程序主要有德国AREVA GmbH的NUDUNA程序[5]、瑞士PSI的NUSS程序[6]、德国GRS的XSUSA程序[7]、比利时SCK·CEN的SANDY程序[8]以及荷兰NRG提出的全蒙特卡罗方法(TMC)[9]。国内西安交通大学基于广义微扰理论开发了燃耗计算灵敏度与不确定度分析程序NECP-SUNDEW[10-11],能够同时量化反应截面、半衰期和裂变产额不确定度在燃耗计算中的传递,清华大学RMC程序实现了随机抽样方法的燃耗不确定度分析模块[12],但只能针对反应截面进行不确定度分析。
本文采用随机抽样方法,以乏燃料同位素成分数据库SFCOMPO-2.0[13]中Takahama-3压水堆组件基准题的SF95-4样品为对象,进行反应截面、衰变常量和裂变产额不确定度引起的核素积存量不确定度的量化分析。
本文采用基于抽样方法的不确定度分析方法。其基本原理就是利用核数据的协方差矩阵进行随机抽样,获得输入参数(核数据)的样本集合,采用输入参数样本计算获得对应的响应结果,并统计产生响应的不确定度。
假设由所需抽样的m个核数据构成的向量为X,其相对协方差矩阵为Σr,先对其进行特征值分解:
Σr=A×Er×AT
(1)
式中:A为Σr特征向量组成的正交矩阵;Er为对角元为Σr特征值的对角矩阵。
(2)
假设向量X满足多元正态分布,则可通过如下随机抽样得到扰动因子向量P[14]:
(3)
式中:I为元素全为1的向量;N(0,1)为由m个满足标准正态分布随机变量组成的随机向量。通过扰动因子向量P可得到X的样本空间:
Xs=PTμ
(4)
式中,μ为X的均值向量。
假设计算样本空间Xs的样本数为N,将计算样本空间Xs作为系统输入,计算不同输入参数样本条件下的响应结果Ri(i=1,2,…,N),对响应R进行统计分析,可得到其期望和标准差信息:
(5)
(6)
式中:E(R)为响应R的期望值;V(R)为响应R的标准差,即不确定度。对于确定论程序,式(6)计算得到的不确定度V(R)就是核数据的不确定度造成响应的不确定度σND。而对于蒙特卡罗程序,每次运行的计算结果都包含蒙特卡罗的统计涨落,即Ri+σMC,i,其中σMC,i为第i次运行结果Ri的蒙特卡罗统计标准差,此时来自核数据的不确定度为:
(7)
根据Wilks公式[15-16],所需要的最小样本数N与概率a以及置信度b的关系满足下式:
1-aN-N(1-a)aN-1≥b
(8)
本文选用300个样本,根据式(8),该样本能覆盖变量的概率密度函数a=98%的可能性为b=98.338 7%,或者a=99%的可能性为b=80.235 0%。
本文使用的反应截面协方差数据来自SCALE6.2软件包中的56群协方差数据库56groupcov,共包含456种不同核素[17]。该数据库中的数据既有来自于ENDF/B-Ⅶ.1的协方差数据,也有美国洛斯阿拉莫斯国家实验室、美国布鲁克海文国家实验室和美国橡树岭国家实验室合作产生的近似数据。
不同核素之间的衰变常量是相互独立的,所以衰变常量的协方差矩阵只有对角元为非零,对角元的值是衰变常量的相对不确定度,该数据直接取自于ENDF/B-Ⅷ.0的MF8/MT457。
在评价核数据库中,同一个质量衰变链中各裂变产物的独立裂变产额之间具有比较大的负相关性。如果直接使用裂变产物的独立产额的不确定度,忽略各裂变产物之间的负相关性,会高估不确定度。本文采用Devillers[18]给出的最小二乘法计算式计算独立产额协方差矩阵,矩阵的对角线元素μii和非对角线元素μij分别为:
(9)
(10)
式中:σi为第i个裂变产物的独立裂变产额的标准差;σ为质量产额的标准差。
裂变产额的不确定度数据来自ENDF/B-Ⅷ.0数据库,并利用式(9)、(10)处理得到裂变产额协方差矩阵,裂变产额协方差数据包括热中子诱发235U和239Pu裂变。图1示出了热中子诱发235U裂变的裂变产额协方差矩阵,核素索引按ZAM(ZAM=10 000Z+10A+M,Z为核电荷数,A为质量数,M为能态)值从小到大排列。
图1 235U热中子裂变产额协方差矩阵
基于上述不确定度分析方法,开发了燃耗计算不确定度分析程序,流程图如图2所示。
图2 燃耗计算不确定度分析程序流程图
本文采用蒙特卡罗程序MCNP6[19]进行输运/燃耗耦合计算,MCNP6程序使用的核数据库为ACE格式连续能量点截面,而协方差矩阵为多群协方差。为了将多群协方差数据产生的扰动因子应用于连续能量点截面的扰动,假设每个能群的能量点完全相关,即同一个能群内的能点共用一个扰动因子。
由于不同反应道的截面之间存在自洽性,因此对某一反应道的扰动,可能还需要对其他反应道进行相应的调整。如果扰动的是总截面,需要同时将扰动量加于对应的各分反应道中。如果对某分反应道进行了修改,总截面也要进行相应调整。同理,对于平均裂变中子数的扰动,也按类似的自洽原则进行处理。
通过随机抽样产生输入样本空间时,有可能产生负值样本,需要进行特殊处理。对于反应截面,本文对负值样本直接舍弃,重新抽样直到样本数满足要求。对于衰变常量和裂变产额,如果扰动因子P<0,则取P=0;如果P>2,则取P=2。
本文基于乏燃料同位素成分数据库SFCOMPO-2.0中Takahama-3压水堆[20]的SF95-4样品进行不确定度分析。图3示出了Takahama-3压水堆17×17燃料组件NT3G23的几何结构,SF95燃料棒位于17a,半径为0.402 5 cm,燃料棒包壳内径、外径分别为0.411、0.475 cm,SF95-4样品燃耗深度为36.69 GW·d/t。
图3 Takahama-3反应堆燃料组件NT3G23的几何结构
燃耗计算由MCNP6程序完成,分别针对反应截面、衰变常量和裂变产额进行不确定度计算,每种类型数据引起的不确定度计算样本数均为300。
图4 反应截面不确定度引起的核素积存量不确定度
表1列出了衰变常量不确定度引起的核素积存量不确定度的计算结果。可看出,除了151Eu,大多数核素的相对不确定度都很小。这是因为151Eu的直接母核151Sm的衰变常量具有8.9%的标准差。
表1 衰变常量不确定度引起的核素积存量不确定度
图5示出了不考虑相关性时裂变产额不确定度引起的部分裂变产物核的核素积存量不确定度,同时也给出了Fiorito等利用SCALE6.2抽样计算的结果。可看出,两者吻合得非常好。图6比较了考虑相关性和不考虑相关性时裂变产额不确定度引起的裂变产物核的核素积存量不确定度。可看出,考虑裂变产额相关性后相对不确定度明显降低。表2列出了考虑裂变产额相关性时裂变产额不确定度引起的部分锕系核素的积存量不确定度。可以看出,裂变产额不确定度引起锕系核素的积存量相对不确定度较小,均小于0.5%。
图5 裂变产额不确定度引起的裂变产物核素积存量不确定度
图6 考虑相关性和不考虑相关性时裂变产额引起的核素积存量不确定度比较
表2 考虑相关性时裂变产额不确定度引起的锕系核素积存量不确定度
SFCOMPO-2.0数据库给出了SF95-4样品的实验测量值,本文对计算值和实验值的偏差与不确定度进行分析。图7示出了SF95-4样品计算值与实验值的相对偏差及其不确定度,相对偏差及其不确定度的定义分别如式(11)、(12)所示:
(11)
(12)
式中:ΔE为实验值的不确定度;ΔC为反应截面、衰变常量和裂变产额不确定度引起的计算值的不确定度。
从图7可看出,对于某些核素,反应截面、衰变常量和裂变产额不确定度引起的不确定度以及实验不确定度无法覆盖计算值与实验值的偏差。如,125Sb的计算值与实验值相对偏差高达169%,针对这种情况,中国核数据中心利用非破坏性实验结果[23]对125Sb实验值进行了分析,怀疑实验数据有问题并进行了修正,修正后的实验值与计算值相对偏差为22%。同理,其他相对偏差比较大的核素,如106Ru、241Am和242Cm的实验值也需要修正实验值。
图7 SF95-4样品计算值与实验值的相对偏差
本文基于抽样方法研究了应用于燃耗计算的不确定度分析方法,研制了不确定度分析程序。基于不确定度分析程序,对Takahama-3压水堆组件基准题中SF95-4样品进行不确定度计算分析,量化了反应截面、衰变常量和裂变产额的不确定度引起的核素积存量不确定度。计算结果表明,本文研制的燃耗计算不确定度分析程序是可靠的,反应截面的不确定度是锕系核积存量不确定度的主要来源,裂变产额和衰变常量的不确定度对部分裂变产物的核素积存量会引起较大的不确定度。但考虑裂变产额相关性后,裂变产额引起的不确定度显著降低。