基于統(tǒng)計(jì)抽樣的敏感性系數(shù)計(jì)算方法

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

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

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

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

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

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

(1)

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

(2)

對(duì)響應(yīng)量進(jìn)行泰勒展開(kāi)，且在一階近似下可得：

(3)

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

(4)

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

var(R)=SVαST

(5)

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

(6)

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

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

(7)

(8)

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

(9)

(10)

上式可寫為：

(11)

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

(12)

2 測(cè)試和驗(yàn)證

2.1 雙群臨界公式驗(yàn)證

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

(13)

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

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

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

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

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

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

(14)

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

2.2 壓水堆單柵元問(wèn)題研究

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

(15)

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

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

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

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

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

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

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

3 結(jié)論

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