韓豐林，祖鐵軍，吳宏春，曹良志

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

輸入?yún)?shù)、數(shù)學(xué)物理模型近似和數(shù)值離散方法等引入的不確定度會(huì)影響核反應(yīng)堆物理計(jì)算結(jié)果的精度，隨著反應(yīng)堆物理計(jì)算方法的發(fā)展，數(shù)學(xué)物理模型近似及數(shù)值離散方法所引入的不確定度顯著降低，而輸入?yún)?shù)對計(jì)算結(jié)果的影響變得更加重要[1-3]。核數(shù)據(jù)是反應(yīng)堆物理計(jì)算中最基礎(chǔ)的輸入?yún)?shù)，它主要是基于核物理實(shí)驗(yàn)測量、核反應(yīng)模型理論分析、核數(shù)據(jù)評價(jià)建庫及宏觀檢驗(yàn)產(chǎn)生的。由于核物理測量的不確定度、核反應(yīng)理論模型的局限性和近似、核數(shù)據(jù)評價(jià)的主觀性、群常數(shù)加工制作中的近似等因素，核工程應(yīng)用中的核數(shù)據(jù)存在一定的不確定度。

核反應(yīng)堆燃耗計(jì)算得到的原子核密度、有效增殖因數(shù)等響應(yīng)參數(shù)，在核反應(yīng)堆設(shè)計(jì)、安全分析及乏燃料處置中十分重要。分析核數(shù)據(jù)的不確定度在燃耗計(jì)算中的傳遞，量化燃耗計(jì)算宏觀響應(yīng)的不確定度，對提高核反應(yīng)堆的安全性、經(jīng)濟(jì)性及高放射性廢物的妥善處置有著積極意義。另外，目前評價(jià)核數(shù)據(jù)庫中的裂變產(chǎn)額數(shù)據(jù)缺乏足夠的獨(dú)立產(chǎn)額測量數(shù)據(jù)和協(xié)方差數(shù)據(jù)[4]，通過靈敏度分析確定對目標(biāo)參數(shù)計(jì)算精度影響較大的裂變產(chǎn)額數(shù)據(jù)，對核數(shù)據(jù)的改進(jìn)有指導(dǎo)意義[5]。

目前國內(nèi)外基于抽樣方法及微擾方法對燃耗計(jì)算響應(yīng)參數(shù)的靈敏度和不確定度開展了一定研究[6-7]，但多針對核反應(yīng)截面開展，缺乏對裂變產(chǎn)額和半衰期等核數(shù)據(jù)的分析。2013年，國際核數(shù)據(jù)評價(jià)合作工作組(WPEC)把量化裂變產(chǎn)額對響應(yīng)參數(shù)引入的不確定度作為未來工作目標(biāo)之一。對于多輸入?yún)?shù)、少輸出參數(shù)的數(shù)值模擬系統(tǒng)，微擾方法與抽樣方法相比計(jì)算效率較高，且可給出響應(yīng)參數(shù)對各輸入?yún)?shù)的靈敏度系數(shù)。

本文基于廣義微擾理論推導(dǎo)燃耗計(jì)算參數(shù)對于裂變產(chǎn)額和半衰期的相對靈敏度系數(shù)計(jì)算理論模型，并在西安交通大學(xué)自主研發(fā)的程序NECP-SUNDEW[8-9]基礎(chǔ)上實(shí)現(xiàn)有效增殖因數(shù)和原子核密度等宏觀響應(yīng)對裂變產(chǎn)額和半衰期的靈敏度計(jì)算，用直接數(shù)值擾動(dòng)方法進(jìn)行驗(yàn)證。以UAM[10-11]基準(zhǔn)題的TMI-1柵元為對象，基于ENDF/B-Ⅶ.1數(shù)據(jù)庫量化由裂變產(chǎn)額和半衰期引入的響應(yīng)參數(shù)不確定度。

1 理論模型

微擾方法量化不確定度首先需計(jì)算出系統(tǒng)的輸入?yún)?shù)對于輸出響應(yīng)參數(shù)的靈敏度，然后基于靈敏度和輸入?yún)?shù)的協(xié)方差數(shù)據(jù)采用“三明治法則”進(jìn)行不確定度計(jì)算，即可獲取輸出參數(shù)的不確定度[12-13]。不確定度的計(jì)算公式為：

(1)

式中：SD(R)為輸出參數(shù)的相對不確定度；Sσ為響應(yīng)參數(shù)R對于核數(shù)據(jù)σ的相對靈敏度系數(shù)矩陣；COVσ為核數(shù)據(jù)σ的相對協(xié)方差矩陣。

1.1 燃耗靈敏度計(jì)算方法

對于系統(tǒng)的響應(yīng)參數(shù)R，如果它是輸入?yún)?shù)α的函數(shù)，輸入?yún)?shù)α的改變造成的響應(yīng)參數(shù)R的相對變化量被定義為響應(yīng)參數(shù)R關(guān)于輸入?yún)?shù)α的相對靈敏度系數(shù)S，其表達(dá)式為：

(2)

在核反應(yīng)堆物理計(jì)算中，響應(yīng)參數(shù)R可表達(dá)為：

Φ*(r,E,Ω,t))drdEdΩdt

(3)

式中：N為原子核密度；Φ為中子角通量密度；Φ*為共軛中子通量密度；r為位置；E為能量；Ω為角度；t為時(shí)間。

將式(2)關(guān)于核數(shù)據(jù)σ進(jìn)行一階泰勒展開：

(4)

式中：t0為燃耗起始時(shí)刻；tf為燃耗的結(jié)束時(shí)刻；t為時(shí)間。式(4)右邊第1項(xiàng)是靈敏度系數(shù)的直接影響項(xiàng)，是由核數(shù)據(jù)的變化直接引起的響應(yīng)參數(shù)的變化；后3項(xiàng)是間接影響項(xiàng)，是由核數(shù)據(jù)的變化引起的各核素的原子核密度、中子通量密度、共軛中子通量密度的變化，并最終造成響應(yīng)參數(shù)的變化。由式(4)可知，求解靈敏度系數(shù)S需給出中子通量密度Φ、共軛中子通量密度Φ*及原子核密度N對于核數(shù)據(jù)的微分。燃耗計(jì)算的基本方程為：

(5)

Bi,g(r)ψi,g(r,Ω)=0

(6)

(7)

(8)

將式(5)～(8)關(guān)于核數(shù)據(jù)σ求微分，可得到計(jì)算式所需的各微分項(xiàng)。若核數(shù)據(jù)為裂變產(chǎn)額和半衰期，經(jīng)推導(dǎo)可得到相對靈敏度系數(shù)計(jì)算式：

(9)

式中：I為燃耗總步數(shù)；N*為共軛原子核密度，由下式得到：

(10)

式中，MT為燃耗矩陣的轉(zhuǎn)置矩陣。

將燃耗矩陣中與半衰期和單群中子通量相關(guān)的部分分開，可寫成以下形式：

M=Mλ+Mφ

(11)

式中：下標(biāo)λ表示燃耗矩陣中由核素衰變引起的轉(zhuǎn)換；下標(biāo)φ表示燃耗矩陣中由中子引起的轉(zhuǎn)換。

結(jié)合式(11)計(jì)算式(9)中的微分項(xiàng)可得到相對靈敏度系數(shù)。

1.2 裂變產(chǎn)額協(xié)方差計(jì)算方法

獨(dú)立裂變產(chǎn)額為裂變直接產(chǎn)生的裂變產(chǎn)物的產(chǎn)額，具有相同質(zhì)量數(shù)的若干裂變產(chǎn)物組成1條質(zhì)量鏈。評價(jià)核數(shù)據(jù)庫中給出的各獨(dú)立裂變產(chǎn)額處于同一條質(zhì)量鏈時(shí)，相互之間有較大的負(fù)相關(guān)性[14]。直接使用評價(jià)核數(shù)據(jù)庫中給出的各獨(dú)立裂變產(chǎn)額標(biāo)準(zhǔn)差數(shù)據(jù)進(jìn)行不確定度計(jì)算，則忽視了獨(dú)立裂變產(chǎn)額的相關(guān)性，會(huì)明顯高估裂變產(chǎn)額引入的不確定度。本文采用Katakura[15]給出的最小二乘法計(jì)算式產(chǎn)生獨(dú)立裂變產(chǎn)額協(xié)方差矩陣，可考慮同一質(zhì)量鏈中各獨(dú)立裂變產(chǎn)額的相關(guān)性，裂變產(chǎn)額協(xié)方差矩陣的對角線元素μii和非對角線元素μij分別為：

(12)

(13)

式中：Δyi為獨(dú)立裂變產(chǎn)額yi的標(biāo)準(zhǔn)差；ΔY為質(zhì)量產(chǎn)額的標(biāo)準(zhǔn)差；j為同一質(zhì)量鏈中的各裂變產(chǎn)物。

在反應(yīng)堆物理計(jì)算程序的燃耗計(jì)算中，一般不采用包含評價(jià)核數(shù)據(jù)庫中全部核素的精細(xì)燃耗鏈，而是在不影響計(jì)算精度的前提下，采用壓縮燃耗鏈以降低計(jì)算資源的需求量。獨(dú)立裂變產(chǎn)額協(xié)方差矩陣包括評價(jià)核數(shù)據(jù)庫中全部裂變產(chǎn)物，需處理得到與壓縮燃耗數(shù)據(jù)庫相一致的裂變產(chǎn)額協(xié)方差矩陣。壓縮燃耗數(shù)據(jù)庫中的裂變產(chǎn)額由精細(xì)燃耗數(shù)據(jù)庫中的獨(dú)立裂變產(chǎn)額按照衰變分支比累積得到，所以壓縮燃耗數(shù)據(jù)庫的裂變產(chǎn)額協(xié)方差矩陣可表示為：

(14)

(15)

式中：A和B均為壓縮燃耗數(shù)據(jù)庫中的裂變產(chǎn)額；y為獨(dú)立裂變產(chǎn)額；a和b為由獨(dú)立裂變產(chǎn)額到壓縮燃耗數(shù)據(jù)庫中裂變產(chǎn)額的轉(zhuǎn)化關(guān)系。

1.3 程序開發(fā)

2 數(shù)值結(jié)果

本文計(jì)算的問題是世界經(jīng)濟(jì)合作與發(fā)展組織核能機(jī)構(gòu)(OECD/NEA)的不確定度分析基準(zhǔn)題(UAM)中的TMI-1柵元，計(jì)算基于熱態(tài)滿功率條件，燃耗期間的平均功率為33.58 W/gU?；鶞?zhǔn)題的幾何結(jié)構(gòu)如圖1所示，表1列出材料信息。本文輸運(yùn)計(jì)算使用69群的多群截面數(shù)據(jù)庫，燃耗計(jì)算使用基于定量化分析方法[17]制作的適用于壓水堆物理計(jì)算的壓縮燃耗數(shù)據(jù)庫，其中包括195種裂變產(chǎn)物。裂變產(chǎn)額和半衰期的不確定度數(shù)據(jù)來自ENDF/B-Ⅶ.1數(shù)據(jù)庫，并經(jīng)過處理得到壓縮燃耗數(shù)據(jù)庫的裂變產(chǎn)額協(xié)方差矩陣，其中裂變產(chǎn)額協(xié)方差數(shù)據(jù)考慮的裂變核素包括235U和239Pu。圖2示出壓縮燃耗數(shù)據(jù)庫的235U熱中子裂變產(chǎn)額協(xié)方差矩陣。

圖1 TMI-1柵元幾何結(jié)構(gòu)Fig.1 Geometry structure of TMI-1 pin-cell

表1 TMI-1柵元材料信息Table 1 Composition of material of TMI-1 pin-cell

圖2 235U熱中子裂變產(chǎn)額協(xié)方差矩陣Fig.2 Fission yield covariance matrix for 235U thermal neutron

2.1 靈敏度分析

使用直接數(shù)值擾動(dòng)方法對基于廣義微擾理論計(jì)算的燃耗相對靈敏度系數(shù)進(jìn)行了驗(yàn)證。表2、3分別列出TMI-1柵元的無限增殖因數(shù)kinf及135I原子核密度對于裂變產(chǎn)額和半衰期的相對靈敏度系數(shù)，其中核數(shù)據(jù)235U_135I表示235U裂變產(chǎn)生135I的產(chǎn)額，GPT和DP分別表示廣義微擾理論和直接數(shù)值擾動(dòng)方法。結(jié)果表明，廣義微擾理論與直接數(shù)值擾動(dòng)方法計(jì)算得到的相對靈敏度系數(shù)結(jié)果偏差很小，證明了開發(fā)的燃耗靈敏度系數(shù)計(jì)算程序的正確性。由相對靈敏度系數(shù)計(jì)算結(jié)果可知，柵元的kinf對裂變產(chǎn)額和半衰期等核數(shù)據(jù)的相對靈敏度系數(shù)較小；而135I原子核密度受裂變產(chǎn)額和半衰期影響明顯，具有較大的相對靈敏度系數(shù)。

原子核密度對于裂變產(chǎn)額和半衰期的相對靈敏度系數(shù)由具體的燃耗轉(zhuǎn)化關(guān)系決定。以135Xe的原子核密度為例，圖3示出其相對靈敏度系數(shù)，圖4示出燃耗計(jì)算中135Xe的燃耗鏈。135Xe的原子核密度對135I裂變產(chǎn)額及135Xe半衰期的相對靈敏度系數(shù)較大，這是由于135Xe本身的裂變產(chǎn)額較小，其產(chǎn)生主要來自135I的衰變，而135Xe的衰變使其消失。隨燃耗深度的增加，235U裂變產(chǎn)額的相對靈敏度系數(shù)變小，而239Pu裂變產(chǎn)額的相對靈敏度系數(shù)變大，這是由于燃耗過程中235U的消耗和239Pu的累積造成的。

表2 kinf的相對靈敏度系數(shù)Table 2 Relative sensitivity coefficient of kinf

表3 135I原子核密度的相對靈敏度系數(shù)Table 3 Relative sensitivity coefficient of 135I nuclide concentration

圖3 135Xe原子核密度的相對靈敏度系數(shù)Fig.3 Relative sensitivity coefficient of 135Xe nuclide concentration

圖4 135Xe燃耗鏈Fig.4 135Xe burnup chain

2.2 不確定度分析

圖5 裂變產(chǎn)額引入的原子核密度的相對不確定度Fig.5 Relative uncertainty of nuclide concentration from fission yield

基于靈敏度分析的結(jié)果，量化了裂變產(chǎn)額和半衰期數(shù)據(jù)對TMI-1柵元的無限增殖因數(shù)kinf和原子核密度計(jì)算結(jié)果引入的相對不確定度。圖5示出壽期末(燃耗深度為60 GW·d/tU)一些核素的原子核密度由裂變產(chǎn)額引入的相對不確定度。結(jié)果表明，考慮裂變產(chǎn)額相關(guān)性后相對不確定度的計(jì)算結(jié)果顯著降低。對于柵元kinf，半衰期引入的相對不確定度極小可忽略，在壽期末僅為0.001%。裂變產(chǎn)額引入的相對不確定度在燃耗過程中的變化如圖6所示，考慮裂變產(chǎn)額相關(guān)性后相對不確定度同樣顯著降低，在壽期末僅為0.06%。對于原子核密度，表4列出壽期末重要裂變產(chǎn)物的原子核密度由裂變產(chǎn)額(考慮相關(guān)性)和半衰期數(shù)據(jù)引入的相對不確定度。結(jié)果表明：半衰期引入的相對不確定度均較小，但151Eu除外，這是由于151Eu的原子核密度對151Sm半衰期的靈敏度很大，而151Sm半衰期具有8.9%的標(biāo)準(zhǔn)差；而裂變產(chǎn)額引入的相對不確定度相對較大，其中109Ag原子核密度的相對不確定度達(dá)到9.9%，主要是由于對應(yīng)的獨(dú)立裂變產(chǎn)額和質(zhì)量產(chǎn)額數(shù)據(jù)有較大的相對不確定度。對于壽期末重核的原子核密度，裂變產(chǎn)額(考慮相關(guān)性)和半衰期數(shù)據(jù)引入的相對不確定度均在0.1%以下。本文對TMI-1柵元的kinf和原子核密度相對不確定度的計(jì)算結(jié)果與Cabellos等[18]使用統(tǒng)計(jì)學(xué)抽樣方法和精細(xì)燃耗鏈獲得的相對不確定度結(jié)果符合。

圖6 裂變產(chǎn)額引入的kinf的相對不確定度Fig.6 Relative uncertainty of kinf from fission yield

表4 裂變產(chǎn)物原子核密度的相對不確定度Table 4 Relative uncertainty of fission product nuclide concentration

3 小結(jié)

基于廣義微擾理論推導(dǎo)了燃耗計(jì)算響應(yīng)參數(shù)對于裂變產(chǎn)額和半衰期的靈敏度計(jì)算模型，在NECP-SUNDEW的基礎(chǔ)上開發(fā)了裂變產(chǎn)額和半衰期的靈敏度和不確定度分析功能。使用開發(fā)的程序計(jì)算了UAM基準(zhǔn)題TMI-1柵元kinf和原子核密度對各種裂變產(chǎn)額和半衰期數(shù)據(jù)的靈敏度，并用直接數(shù)值擾動(dòng)方法對計(jì)算結(jié)果進(jìn)行了驗(yàn)證。發(fā)現(xiàn)柵元kinf對于裂變產(chǎn)額和半衰期的靈敏度較小，而原子核密度對于裂變產(chǎn)額和半衰期的靈敏度較大，由具體的燃耗轉(zhuǎn)化關(guān)系決定?？紤]裂變產(chǎn)額的相關(guān)性，產(chǎn)生了針對壓縮燃耗數(shù)據(jù)庫的裂變產(chǎn)額協(xié)方差矩陣，并量化了柵元kinf和原子核密度由裂變產(chǎn)額和半衰期引入的相對不確定度。柵元kinf由半衰期引入的相對不確定度可忽略，由裂變產(chǎn)額引入的相對不確定度很小，在壽期末僅為0.06%。裂變產(chǎn)額和半衰期對部分裂變產(chǎn)物的原子核密度會(huì)引入較大的相對不確定度，如109Ag和151Eu。