基于抽樣方法的燃耗計(jì)算核素積存量不確定度分析

吳小飛，吳海成，劉 萍，葛智剛

(中國(guó)原子能科學(xué)研究院 核數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室，北京 102413)

燃耗計(jì)算中的核素積存量等響應(yīng)參數(shù)是核反應(yīng)堆設(shè)計(jì)中非常重要的積分參數(shù)，其準(zhǔn)確度直接影響核裝置設(shè)計(jì)的精度。燃耗計(jì)算響應(yīng)參數(shù)的不確定度主要來(lái)源包括輸入?yún)?shù)、物理模型和數(shù)值求解方法。隨著核裝置物理計(jì)算模型和計(jì)算機(jī)技術(shù)的發(fā)展，物理模型簡(jiǎn)化近似和數(shù)值求解方法引入的不確定度顯著降低，輸入?yún)?shù)引入的不確定度已經(jīng)成為響應(yīng)參數(shù)不確定度的主要來(lái)源[1]。核數(shù)據(jù)是燃耗計(jì)算最基本的輸入?yún)?shù)，量化核數(shù)據(jù)不確定性引起的燃耗計(jì)算響應(yīng)參數(shù)的不確定度，對(duì)提高反應(yīng)堆設(shè)計(jì)和乏燃料后處理/貯存/運(yùn)輸?shù)陌踩浴⒔?jīng)濟(jì)性有著重要作用。

目前國(guó)內(nèi)外開(kāi)展的燃耗計(jì)算不確定度分析方法主要分為微擾方法和隨機(jī)抽樣方法。基于燃耗微擾理論[2]的程序較少，主要有日本北海道大學(xué)的CBZ程序[3]、日本JAEA的MARBLE程序[4]等。基于隨機(jī)抽樣方法的程序主要有德國(guó)AREVA GmbH的NUDUNA程序[5]、瑞士PSI的NUSS程序[6]、德國(guó)GRS的XSUSA程序[7]、比利時(shí)SCK·CEN的SANDY程序[8]以及荷蘭NRG提出的全蒙特卡羅方法(TMC)[9]。國(guó)內(nèi)西安交通大學(xué)基于廣義微擾理論開(kāi)發(fā)了燃耗計(jì)算靈敏度與不確定度分析程序NECP-SUNDEW[10-11]，能夠同時(shí)量化反應(yīng)截面、半衰期和裂變產(chǎn)額不確定度在燃耗計(jì)算中的傳遞，清華大學(xué)RMC程序?qū)崿F(xiàn)了隨機(jī)抽樣方法的燃耗不確定度分析模塊[12]，但只能針對(duì)反應(yīng)截面進(jìn)行不確定度分析。

本文采用隨機(jī)抽樣方法，以乏燃料同位素成分?jǐn)?shù)據(jù)庫(kù)SFCOMPO-2.0[13]中Takahama-3壓水堆組件基準(zhǔn)題的SF95-4樣品為對(duì)象，進(jìn)行反應(yīng)截面、衰變常量和裂變產(chǎn)額不確定度引起的核素積存量不確定度的量化分析。

1 不確定度分析

1.1 抽樣方法基本原理

本文采用基于抽樣方法的不確定度分析方法。其基本原理就是利用核數(shù)據(jù)的協(xié)方差矩陣進(jìn)行隨機(jī)抽樣，獲得輸入?yún)?shù)(核數(shù)據(jù))的樣本集合，采用輸入?yún)?shù)樣本計(jì)算獲得對(duì)應(yīng)的響應(yīng)結(jié)果，并統(tǒng)計(jì)產(chǎn)生響應(yīng)的不確定度。

假設(shè)由所需抽樣的m個(gè)核數(shù)據(jù)構(gòu)成的向量為X，其相對(duì)協(xié)方差矩陣為Σr，先對(duì)其進(jìn)行特征值分解：

Σr=A×Er×AT

(1)

式中：A為Σr特征向量組成的正交矩陣；Er為對(duì)角元為Σr特征值的對(duì)角矩陣。

(2)

假設(shè)向量X滿足多元正態(tài)分布，則可通過(guò)如下隨機(jī)抽樣得到擾動(dòng)因子向量P[14]：

(3)

式中：I為元素全為1的向量；N(0,1)為由m個(gè)滿足標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量組成的隨機(jī)向量。通過(guò)擾動(dòng)因子向量P可得到X的樣本空間：

Xs=PTμ

(4)

式中，μ為X的均值向量。

假設(shè)計(jì)算樣本空間Xs的樣本數(shù)為N，將計(jì)算樣本空間Xs作為系統(tǒng)輸入，計(jì)算不同輸入?yún)?shù)樣本條件下的響應(yīng)結(jié)果Ri(i=1,2,…,N)，對(duì)響應(yīng)R進(jìn)行統(tǒng)計(jì)分析，可得到其期望和標(biāo)準(zhǔn)差信息：

(5)

(6)

式中：E(R)為響應(yīng)R的期望值；V(R)為響應(yīng)R的標(biāo)準(zhǔn)差，即不確定度。對(duì)于確定論程序，式(6)計(jì)算得到的不確定度V(R)就是核數(shù)據(jù)的不確定度造成響應(yīng)的不確定度σND。而對(duì)于蒙特卡羅程序，每次運(yùn)行的計(jì)算結(jié)果都包含蒙特卡羅的統(tǒng)計(jì)漲落，即Ri+σMC,i，其中σMC,i為第i次運(yùn)行結(jié)果Ri的蒙特卡羅統(tǒng)計(jì)標(biāo)準(zhǔn)差，此時(shí)來(lái)自核數(shù)據(jù)的不確定度為：

(7)

根據(jù)Wilks公式[15-16]，所需要的最小樣本數(shù)N與概率a以及置信度b的關(guān)系滿足下式：

1-aN-N(1-a)aN-1≥b

(8)

本文選用300個(gè)樣本，根據(jù)式(8)，該樣本能覆蓋變量的概率密度函數(shù)a=98%的可能性為b=98.338 7%，或者a=99%的可能性為b=80.235 0%。

1.2 協(xié)方差數(shù)據(jù)

本文使用的反應(yīng)截面協(xié)方差數(shù)據(jù)來(lái)自SCALE6.2軟件包中的56群協(xié)方差數(shù)據(jù)庫(kù)56groupcov，共包含456種不同核素[17]。該數(shù)據(jù)庫(kù)中的數(shù)據(jù)既有來(lái)自于ENDF/B-Ⅶ.1的協(xié)方差數(shù)據(jù)，也有美國(guó)洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室、美國(guó)布魯克海文國(guó)家實(shí)驗(yàn)室和美國(guó)橡樹嶺國(guó)家實(shí)驗(yàn)室合作產(chǎn)生的近似數(shù)據(jù)。

不同核素之間的衰變常量是相互獨(dú)立的，所以衰變常量的協(xié)方差矩陣只有對(duì)角元為非零，對(duì)角元的值是衰變常量的相對(duì)不確定度，該數(shù)據(jù)直接取自于ENDF/B-Ⅷ.0的MF8/MT457。

在評(píng)價(jià)核數(shù)據(jù)庫(kù)中，同一個(gè)質(zhì)量衰變鏈中各裂變產(chǎn)物的獨(dú)立裂變產(chǎn)額之間具有比較大的負(fù)相關(guān)性。如果直接使用裂變產(chǎn)物的獨(dú)立產(chǎn)額的不確定度，忽略各裂變產(chǎn)物之間的負(fù)相關(guān)性，會(huì)高估不確定度。本文采用Devillers[18]給出的最小二乘法計(jì)算式計(jì)算獨(dú)立產(chǎn)額協(xié)方差矩陣，矩陣的對(duì)角線元素μii和非對(duì)角線元素μij分別為：

(9)

(10)

式中：σi為第i個(gè)裂變產(chǎn)物的獨(dú)立裂變產(chǎn)額的標(biāo)準(zhǔn)差；σ為質(zhì)量產(chǎn)額的標(biāo)準(zhǔn)差。

裂變產(chǎn)額的不確定度數(shù)據(jù)來(lái)自ENDF/B-Ⅷ.0數(shù)據(jù)庫(kù)，并利用式(9)、(10)處理得到裂變產(chǎn)額協(xié)方差矩陣，裂變產(chǎn)額協(xié)方差數(shù)據(jù)包括熱中子誘發(fā)235U和239Pu裂變。圖1示出了熱中子誘發(fā)235U裂變的裂變產(chǎn)額協(xié)方差矩陣，核素索引按ZAM(ZAM=10 000Z+10A+M，Z為核電荷數(shù)，A為質(zhì)量數(shù)，M為能態(tài))值從小到大排列。

圖1 235U熱中子裂變產(chǎn)額協(xié)方差矩陣

1.3 計(jì)算程序研制

基于上述不確定度分析方法，開(kāi)發(fā)了燃耗計(jì)算不確定度分析程序，流程圖如圖2所示。

圖2 燃耗計(jì)算不確定度分析程序流程圖

本文采用蒙特卡羅程序MCNP6[19]進(jìn)行輸運(yùn)/燃耗耦合計(jì)算，MCNP6程序使用的核數(shù)據(jù)庫(kù)為ACE格式連續(xù)能量點(diǎn)截面，而協(xié)方差矩陣為多群協(xié)方差。為了將多群協(xié)方差數(shù)據(jù)產(chǎn)生的擾動(dòng)因子應(yīng)用于連續(xù)能量點(diǎn)截面的擾動(dòng)，假設(shè)每個(gè)能群的能量點(diǎn)完全相關(guān)，即同一個(gè)能群內(nèi)的能點(diǎn)共用一個(gè)擾動(dòng)因子。

由于不同反應(yīng)道的截面之間存在自洽性，因此對(duì)某一反應(yīng)道的擾動(dòng)，可能還需要對(duì)其他反應(yīng)道進(jìn)行相應(yīng)的調(diào)整。如果擾動(dòng)的是總截面，需要同時(shí)將擾動(dòng)量加于對(duì)應(yīng)的各分反應(yīng)道中。如果對(duì)某分反應(yīng)道進(jìn)行了修改，總截面也要進(jìn)行相應(yīng)調(diào)整。同理，對(duì)于平均裂變中子數(shù)的擾動(dòng)，也按類似的自洽原則進(jìn)行處理。

通過(guò)隨機(jī)抽樣產(chǎn)生輸入樣本空間時(shí)，有可能產(chǎn)生負(fù)值樣本，需要進(jìn)行特殊處理。對(duì)于反應(yīng)截面，本文對(duì)負(fù)值樣本直接舍棄，重新抽樣直到樣本數(shù)滿足要求。對(duì)于衰變常量和裂變產(chǎn)額，如果擾動(dòng)因子P<0，則取P=0；如果P>2，則取P=2。

2 數(shù)值計(jì)算與分析

2.1 裝置簡(jiǎn)介

本文基于乏燃料同位素成分?jǐn)?shù)據(jù)庫(kù)SFCOMPO-2.0中Takahama-3壓水堆[20]的SF95-4樣品進(jìn)行不確定度分析。圖3示出了Takahama-3壓水堆17×17燃料組件NT3G23的幾何結(jié)構(gòu)，SF95燃料棒位于17a，半徑為0.402 5 cm，燃料棒包殼內(nèi)徑、外徑分別為0.411、0.475 cm，SF95-4樣品燃耗深度為36.69 GW·d/t。

圖3 Takahama-3反應(yīng)堆燃料組件NT3G23的幾何結(jié)構(gòu)

2.2 計(jì)算結(jié)果與分析

燃耗計(jì)算由MCNP6程序完成，分別針對(duì)反應(yīng)截面、衰變常量和裂變產(chǎn)額進(jìn)行不確定度計(jì)算，每種類型數(shù)據(jù)引起的不確定度計(jì)算樣本數(shù)均為300。

圖4 反應(yīng)截面不確定度引起的核素積存量不確定度

表1列出了衰變常量不確定度引起的核素積存量不確定度的計(jì)算結(jié)果。可看出，除了151Eu，大多數(shù)核素的相對(duì)不確定度都很小。這是因?yàn)?51Eu的直接母核151Sm的衰變常量具有8.9%的標(biāo)準(zhǔn)差。

表1 衰變常量不確定度引起的核素積存量不確定度

圖5示出了不考慮相關(guān)性時(shí)裂變產(chǎn)額不確定度引起的部分裂變產(chǎn)物核的核素積存量不確定度，同時(shí)也給出了Fiorito等利用SCALE6.2抽樣計(jì)算的結(jié)果。可看出，兩者吻合得非常好。圖6比較了考慮相關(guān)性和不考慮相關(guān)性時(shí)裂變產(chǎn)額不確定度引起的裂變產(chǎn)物核的核素積存量不確定度。可看出，考慮裂變產(chǎn)額相關(guān)性后相對(duì)不確定度明顯降低。表2列出了考慮裂變產(chǎn)額相關(guān)性時(shí)裂變產(chǎn)額不確定度引起的部分錒系核素的積存量不確定度。可以看出，裂變產(chǎn)額不確定度引起錒系核素的積存量相對(duì)不確定度較小，均小于0.5%。

圖5 裂變產(chǎn)額不確定度引起的裂變產(chǎn)物核素積存量不確定度

圖6 考慮相關(guān)性和不考慮相關(guān)性時(shí)裂變產(chǎn)額引起的核素積存量不確定度比較

表2 考慮相關(guān)性時(shí)裂變產(chǎn)額不確定度引起的錒系核素積存量不確定度

2.3 討論

SFCOMPO-2.0數(shù)據(jù)庫(kù)給出了SF95-4樣品的實(shí)驗(yàn)測(cè)量值，本文對(duì)計(jì)算值和實(shí)驗(yàn)值的偏差與不確定度進(jìn)行分析。圖7示出了SF95-4樣品計(jì)算值與實(shí)驗(yàn)值的相對(duì)偏差及其不確定度，相對(duì)偏差及其不確定度的定義分別如式(11)、(12)所示：

(11)

(12)

式中：ΔE為實(shí)驗(yàn)值的不確定度；ΔC為反應(yīng)截面、衰變常量和裂變產(chǎn)額不確定度引起的計(jì)算值的不確定度。

從圖7可看出，對(duì)于某些核素，反應(yīng)截面、衰變常量和裂變產(chǎn)額不確定度引起的不確定度以及實(shí)驗(yàn)不確定度無(wú)法覆蓋計(jì)算值與實(shí)驗(yàn)值的偏差。如，125Sb的計(jì)算值與實(shí)驗(yàn)值相對(duì)偏差高達(dá)169%，針對(duì)這種情況，中國(guó)核數(shù)據(jù)中心利用非破壞性實(shí)驗(yàn)結(jié)果[23]對(duì)125Sb實(shí)驗(yàn)值進(jìn)行了分析，懷疑實(shí)驗(yàn)數(shù)據(jù)有問(wèn)題并進(jìn)行了修正，修正后的實(shí)驗(yàn)值與計(jì)算值相對(duì)偏差為22%。同理，其他相對(duì)偏差比較大的核素，如106Ru、241Am和242Cm的實(shí)驗(yàn)值也需要修正實(shí)驗(yàn)值。

圖7 SF95-4樣品計(jì)算值與實(shí)驗(yàn)值的相對(duì)偏差

3 總結(jié)

本文基于抽樣方法研究了應(yīng)用于燃耗計(jì)算的不確定度分析方法，研制了不確定度分析程序。基于不確定度分析程序，對(duì)Takahama-3壓水堆組件基準(zhǔn)題中SF95-4樣品進(jìn)行不確定度計(jì)算分析，量化了反應(yīng)截面、衰變常量和裂變產(chǎn)額的不確定度引起的核素積存量不確定度。計(jì)算結(jié)果表明，本文研制的燃耗計(jì)算不確定度分析程序是可靠的，反應(yīng)截面的不確定度是錒系核積存量不確定度的主要來(lái)源，裂變產(chǎn)額和衰變常量的不確定度對(duì)部分裂變產(chǎn)物的核素積存量會(huì)引起較大的不確定度。但考慮裂變產(chǎn)額相關(guān)性后，裂變產(chǎn)額引起的不確定度顯著降低。