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

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

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

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

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

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

1 理論模型

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

(1)

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

1.1 燃耗靈敏度計算方法

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

(2)

在核反應堆物理計算中，響應參數(shù)R可表達為：

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

(3)

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

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

(4)

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

(5)

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

(6)

(7)

(8)

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

(9)

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

(10)

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

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

M=Mλ+Mφ

(11)

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

結合式(11)計算式(9)中的微分項可得到相對靈敏度系數(shù)。

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

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

(12)

(13)

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

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

(14)

(15)

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

1.3 程序開發(fā)

2 數(shù)值結果

本文計算的問題是世界經(jīng)濟合作與發(fā)展組織核能機構(OECD/NEA)的不確定度分析基準題(UAM)中的TMI-1柵元，計算基于熱態(tài)滿功率條件，燃耗期間的平均功率為33.58 W/gU?；鶞暑}的幾何結構如圖1所示，表1列出材料信息。本文輸運計算使用69群的多群截面數(shù)據(jù)庫，燃耗計算使用基于定量化分析方法[17]制作的適用于壓水堆物理計算的壓縮燃耗數(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柵元幾何結構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ù)值擾動方法對基于廣義微擾理論計算的燃耗相對靈敏度系數(shù)進行了驗證。表2、3分別列出TMI-1柵元的無限增殖因數(shù)kinf及135I原子核密度對于裂變產(chǎn)額和半衰期的相對靈敏度系數(shù)，其中核數(shù)據(jù)235U_135I表示235U裂變產(chǎn)生135I的產(chǎn)額，GPT和DP分別表示廣義微擾理論和直接數(shù)值擾動方法。結果表明，廣義微擾理論與直接數(shù)值擾動方法計算得到的相對靈敏度系數(shù)結果偏差很小，證明了開發(fā)的燃耗靈敏度系數(shù)計算程序的正確性。由相對靈敏度系數(shù)計算結果可知，柵元的kinf對裂變產(chǎn)額和半衰期等核數(shù)據(jù)的相對靈敏度系數(shù)較?。欢?35I原子核密度受裂變產(chǎn)額和半衰期影響明顯，具有較大的相對靈敏度系數(shù)。

原子核密度對于裂變產(chǎn)額和半衰期的相對靈敏度系數(shù)由具體的燃耗轉(zhuǎn)化關系決定。以135Xe的原子核密度為例，圖3示出其相對靈敏度系數(shù)，圖4示出燃耗計算中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

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

圖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 小結

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