馬續(xù)波，唐 輝，楊 樂，朱帥濤，陳義學(xué)

(華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206)

不確定性和敏感性分析在各領(lǐng)域都有重要地位[1-3]，如對物理模型輸入?yún)?shù)的敏感性分析可知哪個參數(shù)的重要性更大，確定出問題的主要因素，以方便構(gòu)造更好的模型。在核工程領(lǐng)域，為保證核電廠的絕對安全，以往的計算模型通常采用留有很大余量的保守模型進行設(shè)計，而進行反應(yīng)堆工程設(shè)計中的關(guān)鍵參數(shù)的不確定性和敏感性分析方法可較好理解核電廠設(shè)計的保守性以及改進核電廠的經(jīng)濟性。以往的不確定性和敏感性分析方法通常有兩種，分別是確定論方法和統(tǒng)計抽樣方法，但基于確定論方法的不確定性分析方法通常僅能考慮低階的不確定性(一般能到二階，高階較難考慮)，并且對分析的響應(yīng)量有一定要求。所以發(fā)展了廣義微擾理論的敏感性計算方法，以可更好考慮響應(yīng)量的擴展量，如控制棒價值、堆芯功率等。而基于統(tǒng)計抽樣方法的不確定性分析方法，其優(yōu)點是對響應(yīng)量無要求，且易計算響應(yīng)量的不確定度，其缺點是較難計算響應(yīng)量輸入?yún)?shù)的敏感性系數(shù)[4-5]以及耗費的計時較長。伴隨計算機計算能力的發(fā)展，計算時間較長問題也在逐步得到改善。本文推導(dǎo)基于統(tǒng)計抽樣方法計算響應(yīng)量對輸入?yún)?shù)的敏感性系數(shù)的理論公式，然后利用算例進行驗證分析，驗證理論方法的可行性。

1 敏感性系數(shù)計算理論

假設(shè)響應(yīng)量R為N個輸入變量的函數(shù)：

R=R(α1,…,αN)

(1)

而N個輸入變量的真值是由變量的最佳估計值和相應(yīng)的不確定度組成的：

(2)

對響應(yīng)量進行泰勒展開，且在一階近似下可得：

(3)

通常定義輸入變量αi相對于響應(yīng)量R的敏感性系數(shù)為：

(4)

如果計算得到了響應(yīng)量R對每個參數(shù)的敏感性系數(shù)，則可根據(jù)誤差傳遞公式(式(5))計算響應(yīng)量R的協(xié)方差，進而利用式(6)求得響應(yīng)量R的標準偏差ΔR，即為R的不確定度[3]。

var(R)=SVαST

(5)

ΔR=[var(R)]1/2

(6)

其中：S為敏感性系數(shù)向量，是N個參數(shù)組成的行向量；ST為向量S的轉(zhuǎn)置向量；Vα為N×N的輸入變量的協(xié)方差矩陣。

通過上面的分析可知：求ΔR需先求得響應(yīng)量R對每個輸入?yún)?shù)αi的敏感性系數(shù)矩陣S，根據(jù)輸入變量之間的協(xié)方差矩陣可求得ΔR。而基于統(tǒng)計抽樣方法的不確定性分析方法的思路與上面有所不同，其基本思想為：先根據(jù)輸入變量α中每個參數(shù)的分布規(guī)律以及每個參數(shù)之間的協(xié)方差矩陣Vα進行抽樣，然后利用式(1)可得到響應(yīng)量R的樣本，對樣本進行統(tǒng)計分析即可得到ΔR?；诮y(tǒng)計抽樣方法的優(yōu)點在于：1) 不僅可考慮輸入?yún)?shù)對R的一階擾動，所有階的擾動均可考慮；2) 對響應(yīng)量無限制，方法的通用性更好。但其缺點也很明顯，即不易求出輸入變量αi相對于響應(yīng)量R的敏感性系數(shù)。在核工程中，由于離散能群的不等寬效應(yīng)，可能產(chǎn)生因截面所屬能量寬度較大引起的較大敏感性系數(shù)[6]，為此，本文引入平均勒能量相對敏感性系數(shù)：

(7)

(8)

式中：i為樣本編號；j為輸入變量即能群截面的編號。上式可變?yōu)椋?/p>

(9)

(10)

上式可寫為：

(11)

利用上式，可求得響應(yīng)量對每個能群截面的敏感性系數(shù)：

(12)

2 測試和驗證

2.1 雙群臨界公式驗證

為驗證基于統(tǒng)計抽樣的敏感性系數(shù)計算方法的可行性，本文提出了利用裸堆雙群近似下的臨界公式進行驗證。裸堆雙群近似的臨界公式如式(13)所示。

(13)

其中：(υΣf)1、(υΣf)2分別為第1和第2能群的產(chǎn)生截面；D1、D2分別為第1和第2能群的擴散系數(shù)；Σ1→2為第1能群到第2能群的散射截面；Σa,1、Σa,2分別為第1、2能群的吸收截面；Σr=Σa,1+Σ1→2，為輸運截面；B2為幾何曲率常數(shù)，假定幾何曲率常數(shù)B2=2.9×10-4。

選取了式中的7個參數(shù)進行研究，7個參數(shù)的編號及初始值列于表1。假定每個參數(shù)的相對誤差為1%，且服從正態(tài)分布。假定7個參數(shù)之間是相互獨立的，不考慮參數(shù)之間的相關(guān)性情況下，表2列出了采用不同的樣本總數(shù)的敏感性系數(shù)的結(jié)果，表2中的基準值是指直接利用式(4)的計算結(jié)果。由表2的結(jié)果可見：基于抽樣方法的敏感性系數(shù)計算方法是正確的，且對于敏感性系數(shù)較大的參數(shù)，隨著抽樣樣本數(shù)的增加，計算結(jié)果與基準值吻合越好，但對于敏感性系數(shù)很小的參數(shù)，如D2和Σr，由于統(tǒng)計性原因?qū)е抡`差較大。圖1a示出了敏感性系數(shù)的相對誤差隨樣本總數(shù)的變化。由圖1a可見，樣本數(shù)達到2 000時，除了兩個敏感性系數(shù)較小的參數(shù)，其他各主要的誤差基本可控制在5%以內(nèi)。這也說明了該方法的可行性。

表1 裸堆雙群近似的臨界公式參數(shù)編號及初始值Table 1 Parameter number and initial value of two groups multiplication factor of bare homogeneous ractor

基于抽樣的敏感性系數(shù)計算方法的1個優(yōu)點是可方便考慮研究參數(shù)之間的相關(guān)性。本文研究了前3個參數(shù)的相關(guān)性對敏感性系數(shù)計算的影響。假定前3個參數(shù)的相關(guān)性系數(shù)矩陣為Σ，有：

表2 不同抽樣樣本數(shù)對結(jié)果影響Table 2 Result effects with respect to sampling number

a——敏感性系數(shù)的相對誤差隨抽樣樣本數(shù)的變化關(guān)系；b——(υΣf)1、(υΣf)2和Σ1→2之間的相關(guān)性系數(shù)對敏感性系數(shù)計算的影響圖1 抽樣樣本數(shù)和參數(shù)之間相關(guān)性對敏感性系數(shù)的影響Fig.1 Sensitivity coefficients variety with total sampling sizes and correlation between parameters

(14)

其中，a為兩個參數(shù)的相關(guān)性的大小。圖1b示出了設(shè)定的相關(guān)性系數(shù)與敏感性系數(shù)的關(guān)系(樣本數(shù)為10 000)，由圖1b可見，在相關(guān)性系數(shù)較小的情況下，相關(guān)性系數(shù)對3個參數(shù)的敏感性系數(shù)的計算結(jié)果影響不大，但在相關(guān)性系數(shù)較大的情況下，有可能導(dǎo)致相關(guān)性系數(shù)出現(xiàn)跳躍現(xiàn)象，且在后面的壓水堆單柵元的敏感性分析中也出現(xiàn)了類似現(xiàn)象。

2.2 壓水堆單柵元問題研究

核工程設(shè)計中，通過對復(fù)雜系統(tǒng)的不確定性分析可減小系統(tǒng)的近似程度，從而提高核設(shè)計的經(jīng)濟性。通過對系統(tǒng)的敏感性系數(shù)分析，可進一步明確輸入變量的重要性，進而為提高輸入?yún)?shù)精度指明方向。近年來，國內(nèi)外對此也進行了大量研究[7-17]。為研究核工程計算過程中的不確定性，OECD/NEA專門提出了基準題，國際上所有研究單位可利用此基準題對比校核各自程序的適用性[9]。本文采用基準題中的三哩島(TMI)燃料柵元為研究對象，研究利用統(tǒng)計抽樣方法計算235U裂變截面的敏感性系數(shù)。為解決式中的協(xié)方差求逆問題，采用兩種方法解決此問題。方法1是去掉協(xié)方差矩陣中的相關(guān)項，只保留主對角線部分進行求逆；方法2是假定每個能群的截面的誤差均為0.1%或1%的協(xié)方差矩陣。

(15)

式中，V′σ,σ為簡化后的協(xié)方差矩陣。

235U裂變截面的協(xié)方差矩陣采用69群能群結(jié)構(gòu)，基于ENDF/B-Ⅶ.1核評價數(shù)據(jù)庫，經(jīng)NJOY程序加工而成[18]。235U裂變截面的相關(guān)性系數(shù)矩陣如圖2所示。由圖2可知，部分能群之間存在很強的相關(guān)性。

為驗證不同計算方法的正確性，開發(fā)不確定性與敏感性分析程序SUACL，該程序目前主要具備壓水堆開展截面的不確定性和敏感性系數(shù)分析功能。計算中對TMI的輸運計算采用了加拿大蒙特利爾大學(xué)開發(fā)的反應(yīng)堆計算模擬軟件DRAGON4.0。采用方法1或方法2在計算敏感性系數(shù)時均采用式(15)，不再采用式(12)，因為產(chǎn)生樣本時的協(xié)方差矩陣已用V′σ,σ替換了Vσ,σ，由于V′σ,σ很易求出其逆矩陣，因此可很方便求出。圖3示出了方法1計算結(jié)果與直接微擾計算結(jié)果對比，圖中的Full matrix表示采用圖2給出的協(xié)方差矩陣，抽樣樣本為500時計算得到的相對敏感性系數(shù)；Diagonal matrix表示對圖2給出的協(xié)方差矩陣去掉了相關(guān)性系數(shù)，僅保留了主對角線部分，即方法1，計算得到的相對敏感性系數(shù)；Direct perturbation表示直接利用式(4)計算得到的結(jié)果。由圖3可見，采用方法1，即主對角元矩陣，而不考慮相關(guān)性，計算結(jié)果與直接微擾的結(jié)果吻合很好。這種結(jié)果很易理解，因為利用式(4)的直接微擾方法中并未考慮能群之間的相關(guān)性。而考慮能群之間的相關(guān)性系數(shù)后會出現(xiàn)震蕩，通過之前的雙群近似下的臨界公式分析可知，震蕩主要是由于考慮了能群之間的相關(guān)性導(dǎo)致。

圖2 235U 69群裂變截面的相關(guān)性系數(shù)矩陣Fig.2 Correlation coefficient matrix of 69 groups fission cross section of 235U

圖3 利用方法1計算的235U裂變截面的相對敏感性系數(shù)Fig.3 Relative sensitivity coefficient of 235U fission cross section with method 1

方法1中采用主對角線元素就可得到較好的相對敏感性系數(shù)，由此推斷：如果每個能群均采用相同的擾動量，如均為1%或0.1%，那么也應(yīng)能計算出較好的敏感性系數(shù)，為此采用方法2進行驗證。圖4示出了方法2中各能群截面的微擾量為1%和0.1%的計算結(jié)果與直接微擾法的比較，樣本總數(shù)為500。由圖4可見，在本例子中，微擾量1%和0.1%均可計算出235U裂變截面的相對敏感性系數(shù)，但微擾量為0.1%的計算結(jié)果較1%的微擾量的計算結(jié)果要好，即與直接微擾法相比，誤差更好。

圖4 利用方法2計算的235U裂變截面的相對敏感性系數(shù)Fig.4 Relative sensitivity cofficient of 235U fission cross section with method 2

根據(jù)敏感性系數(shù)可利用式(5)計算得到對應(yīng)響應(yīng)量的不確定度。表3列出了采用不同的敏感性系數(shù)計算方法下235U裂變截面對柵元kinf的相對不確定度及其誤差。由表3可見，微擾量為0.1%時較微擾量為1%情況好，此結(jié)果也可由圖4看出，采用Full matrix的情況，雖然敏感性系數(shù)跳躍性很強，但很多能群正負抵消，導(dǎo)致最后對柵元的kinf的相對不確定度的誤差僅1.36%。

表3 不同敏感性系數(shù)計算方法下kinf的不確定度Table 3 Uncertainty of kinf with different sensitivity cofficient calculation methods

3 結(jié)論

本文推導(dǎo)了基于抽樣方法進行敏感性系數(shù)計算的理論公式，首先利用簡單的例子進行了驗證分析，驗證了方法的可行性，同時也發(fā)現(xiàn)了該方法在計算敏感性系數(shù)較小的參數(shù)時結(jié)果可能會有較大誤差。針對核工程中的單柵元基準題問題，利用此方法的難點為協(xié)方差矩陣的求逆，本文提出了兩種替代的方法進行敏感性系數(shù)計算，從而避免了協(xié)方差矩陣求逆困難的問題，進而對所提出的方法進行了驗證，驗證了方法的可行性。本文還發(fā)現(xiàn)微擾量的不同會導(dǎo)致計算的敏感性系數(shù)的精度不同，微擾量的選取方法是下一步要研究的問題。