針對裂變產額和半衰期的燃耗計算靈敏度和不確定度分析方法

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

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

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

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

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

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

1 理論模型

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

(1)

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

1.1 燃耗靈敏度計算方法

對于系統的響應參數R，如果它是輸入參數α的函數，輸入參數α的改變造成的響應參數R的相對變化量被定義為響應參數R關于輸入參數α的相對靈敏度系數S，其表達式為：

(2)

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

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

(3)

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

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

(4)

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

(5)

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

(6)

(7)

(8)

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

(9)

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

(10)

式中，MT為燃耗矩陣的轉置矩陣。

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

M=Mλ+Mφ

(11)

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

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

1.2 裂變產額協方差計算方法

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

(12)

(13)

式中：Δyi為獨立裂變產額yi的標準差；ΔY為質量產額的標準差；j為同一質量鏈中的各裂變產物。

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

(14)

(15)

式中：A和B均為壓縮燃耗數據庫中的裂變產額；y為獨立裂變產額；a和b為由獨立裂變產額到壓縮燃耗數據庫中裂變產額的轉化關系。

1.3 程序開發

2 數值結果

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

圖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熱中子裂變產額協方差矩陣Fig.2 Fission yield covariance matrix for 235U thermal neutron

2.1 靈敏度分析

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

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

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

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

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

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

2.2 不確定度分析

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

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

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

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

3 小結

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