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

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

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

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

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

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

1 不確定度分析

1.1 抽樣方法基本原理

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

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

Σr=A×Er×AT

(1)

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

(2)

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

(3)

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

Xs=PTμ

(4)

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

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

(5)

(6)

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

(7)

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

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

(8)

本文選用300個樣本，根據(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ù)來自SCALE6.2軟件包中的56群協(xié)方差數(shù)據(jù)庫56groupcov，共包含456種不同核素[17]。該數(shù)據(jù)庫中的數(shù)據(jù)既有來自于ENDF/B-Ⅶ.1的協(xié)方差數(shù)據(jù)，也有美國洛斯阿拉莫斯國家實驗室、美國布魯克海文國家實驗室和美國橡樹嶺國家實驗室合作產(chǎn)生的近似數(shù)據(jù)。

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

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

(9)

(10)

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

裂變產(chǎn)額的不確定度數(shù)據(jù)來自ENDF/B-Ⅷ.0數(shù)據(jù)庫，并利用式(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 計算程序研制

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

圖2 燃耗計算不確定度分析程序流程圖

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

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

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

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

2.1 裝置簡介

本文基于乏燃料同位素成分?jǐn)?shù)據(jù)庫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 計算結(jié)果與分析

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

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

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

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

圖5示出了不考慮相關(guān)性時裂變產(chǎn)額不確定度引起的部分裂變產(chǎn)物核的核素積存量不確定度，同時也給出了Fiorito等利用SCALE6.2抽樣計算的結(jié)果?？煽闯觯瑑烧呶呛系梅浅：?。圖6比較了考慮相關(guān)性和不考慮相關(guān)性時裂變產(chǎn)額不確定度引起的裂變產(chǎn)物核的核素積存量不確定度?？煽闯?，考慮裂變產(chǎn)額相關(guān)性后相對不確定度明顯降低。表2列出了考慮裂變產(chǎn)額相關(guān)性時裂變產(chǎn)額不確定度引起的部分錒系核素的積存量不確定度?？梢钥闯?，裂變產(chǎn)額不確定度引起錒系核素的積存量相對不確定度較小，均小于0.5%。

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

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

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

2.3 討論

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

(11)

(12)

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

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

圖7 SF95-4樣品計算值與實驗值的相對偏差

3 總結(jié)

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