核截面不確定性引起的keff不確定度的直接蒙特卡羅方法計(jì)算

王新哲，喻 宏，王文明，胡 赟，楊曉燕

(1.中國(guó)原子能科學(xué)研究院 快堆研究設(shè)計(jì)所，北京 102413；2.中國(guó)原子能科學(xué)研究院 核物理研究所，北京 102413)

核裝置的有效增殖因數(shù)keff在核裝置設(shè)計(jì)中起至關(guān)重要的作用，其準(zhǔn)確度直接影響核裝置設(shè)計(jì)的精度，keff計(jì)算的不確定度主要源于以下3個(gè)方面：1) 數(shù)據(jù)的不確定性，包括核數(shù)據(jù)、裝置尺寸以及裝載量等數(shù)據(jù)的不確定度造成的不確定性；2) 計(jì)算模型建立中的簡(jiǎn)化、近似造成的不確定性；3) 計(jì)算方法與計(jì)算過(guò)程中的簡(jiǎn)化、近似造成的不確定性。隨著計(jì)算機(jī)硬件及計(jì)算方法的不斷發(fā)展，模型建立及計(jì)算方法導(dǎo)致keff的不確定度越來(lái)越小，影響核裝置keff不確定度的重要因素在于數(shù)據(jù)的不確定性，這其中又以核截面數(shù)據(jù)的不確定度為主。

現(xiàn)階段確定核截面數(shù)據(jù)不確定性對(duì)核裝置keff計(jì)算不確定度的影響的方法主要為靈敏度-不確定度分析方法(S-U方法)[1]，該方法的主要思想是：首先通過(guò)直接微擾或共軛方法(廣義微擾)確定核裝置keff對(duì)各核素各截面數(shù)據(jù)的靈敏度，然后再與協(xié)方差矩陣進(jìn)行數(shù)學(xué)處理，即可得到核截面數(shù)據(jù)不確定性造成的裝置keff計(jì)算不確定度。該方法對(duì)于一、二維情況，具有計(jì)算速度快并可計(jì)算各種核素各截面靈敏度與不確定度的優(yōu)點(diǎn)，同時(shí)，能進(jìn)行不確定度來(lái)源分析，為改善核數(shù)據(jù)指明方向。但對(duì)于三維情況，該方法需用輸運(yùn)或蒙特卡羅方法計(jì)算通量的價(jià)值函數(shù)，在現(xiàn)階段仍存在計(jì)算時(shí)間長(zhǎng)、計(jì)算精度差、計(jì)算較繁瑣等問(wèn)題。本文利用直接蒙特卡羅方法計(jì)算核截面數(shù)據(jù)不確定性對(duì)核裝置keff計(jì)算不確定度的影響，以能在現(xiàn)有計(jì)算平臺(tái)上簡(jiǎn)單、有效地給出keff計(jì)算不確定度結(jié)果。

1 理論分析

蒙特卡羅方法的主要思想是當(dāng)所求問(wèn)題的解是某事件的概率，或某隨機(jī)變量的數(shù)學(xué)期望，或與概率、數(shù)學(xué)期望有關(guān)的量時(shí)，通過(guò)某種試驗(yàn)的方法，得出該事件發(fā)生的頻率，或該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值，通過(guò)它得到問(wèn)題的解[2]。

直接蒙特卡羅不確定分析方法在反應(yīng)堆物理領(lǐng)域，由于計(jì)算條件及各種核數(shù)據(jù)不確定度之間的強(qiáng)相關(guān)性等原因，尚未采用。

核裝置keff計(jì)算中使用的核截面數(shù)據(jù)，主要由實(shí)驗(yàn)測(cè)量得到及對(duì)測(cè)量數(shù)據(jù)進(jìn)行數(shù)學(xué)加工得到。這兩種過(guò)程中，均會(huì)產(chǎn)生一定的不確定度。同時(shí)，由于實(shí)驗(yàn)測(cè)量中的系統(tǒng)誤差等因素，各核數(shù)據(jù)的不確定度相互關(guān)聯(lián)。在評(píng)價(jià)核數(shù)據(jù)庫(kù)中，用協(xié)方差數(shù)據(jù)的形式來(lái)描述核數(shù)據(jù)的不確定度及其相互關(guān)聯(lián)程度[3]。

直接蒙特卡羅方法求解核截面數(shù)據(jù)不確定性對(duì)核裝置keff影響的基本思路為：1) 產(chǎn)生1組符合所給截面及相應(yīng)協(xié)方差數(shù)據(jù)的隨機(jī)截面數(shù)據(jù)；2) 運(yùn)用該組核截面數(shù)據(jù)通過(guò)堆芯計(jì)算程序進(jìn)行keff計(jì)算，并記錄keff；3) 重復(fù)上述兩步，直到達(dá)到要求的模擬次數(shù)；4) 對(duì)keff進(jìn)行統(tǒng)計(jì)，得出該組數(shù)據(jù)的期望值及標(biāo)準(zhǔn)差，所得結(jié)果即為keff因核截面數(shù)據(jù)不確定性導(dǎo)致的不確定度。計(jì)算流程如圖1所示。

圖1 直接蒙特卡羅方法流程圖

1.1 隨機(jī)截面數(shù)據(jù)產(chǎn)生

從計(jì)算流程可知，隨機(jī)截面數(shù)據(jù)的產(chǎn)生是直接蒙特卡羅方法中最關(guān)鍵的一步。本文計(jì)算過(guò)程中假定截面數(shù)據(jù)的分布符合正態(tài)分布，這就要求產(chǎn)生的截面數(shù)據(jù)既要各自滿(mǎn)足正態(tài)分布標(biāo)準(zhǔn)，同時(shí)又要符合協(xié)方差數(shù)據(jù)要求。由于數(shù)據(jù)限制，本文僅使用各反應(yīng)道自身群截面協(xié)方差數(shù)據(jù)，但只要能獲得相關(guān)數(shù)據(jù)，即可很容易地?cái)U(kuò)展到不同反應(yīng)道之間的計(jì)算。

產(chǎn)生n群隨機(jī)截面數(shù)據(jù)，首先要求產(chǎn)生n維相互無(wú)關(guān)符合標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)向量，該隨機(jī)向量可通過(guò)Box-Muller方法[4]由式(1)兩兩產(chǎn)生：

(1)

其中：U1、U2為在[0，1]上符合均勻分布的隨機(jī)數(shù)；X1、X2為符合(0，1)分布的正態(tài)相互無(wú)關(guān)隨機(jī)數(shù)，通過(guò)上述方法，即可產(chǎn)生n維相互無(wú)關(guān)正態(tài)向量X，X～(0，I)(X符合(0，I)分布，其中，I為n維單位矩陣)。

由正態(tài)變量的線性變換不變性，即n維隨機(jī)向量X～(a，B)，若B為n維正定矩陣，C為m×n階矩陣，且行向量線性無(wú)關(guān)，則m維隨機(jī)向量Y=CX～(Ca,CBCT)，其中，CT表示矩陣C的轉(zhuǎn)置[5]。利用該條性質(zhì)，可將X～(0，I)變換成符合要求的隨機(jī)向量Y～(sig，C)，sig為核數(shù)據(jù)中的截面向量。

由于假設(shè)核截面數(shù)據(jù)符合n維隨機(jī)正態(tài)分布，故其協(xié)方差陣C對(duì)稱(chēng)正定[6]，對(duì)于對(duì)稱(chēng)正定陣，可通過(guò)Cholesky分解得到一個(gè)下三角與一個(gè)上三角且互為轉(zhuǎn)置的兩個(gè)矩陣的乘積[7]，如式(2)所示，其中，D為上三角矩陣。

C=DTD

(2)

則對(duì)X～(0，I)，有：

DX～(DT0,DTID)～(0,C)

(3)

對(duì)DX向量進(jìn)行平移，即得：

Y=sig+DX～(sig,C)

(4)

其中，sig為未經(jīng)抽樣原始截面數(shù)據(jù)。通過(guò)上述過(guò)程，即可產(chǎn)生符合Y～(sig，C)的隨機(jī)截面數(shù)據(jù)。

1.2 計(jì)算結(jié)果分析

本文采用的堆芯計(jì)算程序?yàn)榛谟邢薏罘址椒ǖ腃ITATION程序[8]和基于粗網(wǎng)格節(jié)塊擴(kuò)散方法的NAS程序[9]，其基本計(jì)算原理均源自多群擴(kuò)散方程，故只需討論擴(kuò)散方程中計(jì)算結(jié)果的分布即可。算符形式多群擴(kuò)散方程為：

(5)

其中：M為消失項(xiàng)算符；F為產(chǎn)生項(xiàng)算符；Σf、ΣR分別為宏觀裂變截面和宏觀移出截面；χ為裂變中子譜；ν為有效裂變中子數(shù)。對(duì)于核截面數(shù)據(jù)引起的不確定度分析，可假定除截面外其他各量均為常量，同時(shí)截面向量符合正態(tài)分布。對(duì)于算符F，其中每項(xiàng)均符合正態(tài)分布，則由正態(tài)向量運(yùn)算法則可知，F(xiàn)Φ為正態(tài)隨機(jī)向量。對(duì)于算符M，除泄漏項(xiàng)外各項(xiàng)皆符合正態(tài)分布，泄漏項(xiàng)經(jīng)過(guò)差分后變成擴(kuò)散系數(shù)D的線性運(yùn)算，因此只需討論D的分布即可。

按照擴(kuò)散理論，D的分量Dg為：

(6)

(7)

Dg的分布為倒數(shù)正態(tài)分布，當(dāng)μ≥σ時(shí)，Dg近似符合正態(tài)分布。在核數(shù)據(jù)尤其是少群核數(shù)據(jù)中，該假設(shè)成立，則：

(8)

其中，X、Y均為多維正態(tài)列向量。此時(shí)，式(4)可改寫(xiě)為：

kX=Y

(9)

由擴(kuò)散方程物理意義可知，X、Y中各元素均大于0，故一定存在各項(xiàng)皆為正數(shù)的行向量x-1，使得：

x-1X=1

(10)

將式(9)等號(hào)兩邊均乘x-1，即可得：

k=x-1Y

(11)

則由正態(tài)變量的線性變換不變性，可得：

(12)

因此，可證得k近似符合正態(tài)分布。將程序所得的結(jié)果進(jìn)行正態(tài)擬合，可得到數(shù)學(xué)期望和標(biāo)準(zhǔn)差，所得標(biāo)準(zhǔn)差數(shù)據(jù)即為核截面數(shù)據(jù)不確定性導(dǎo)致的核裝置keff不確定度。

2 程序開(kāi)發(fā)

(13)

3 算法驗(yàn)證及結(jié)果分析

本文分別采用小型核裝置Jezebel-239Pu和中國(guó)實(shí)驗(yàn)快堆(CEFR)首爐裝料作為算法驗(yàn)證模型。

3.1 Jezebel-239Pu裝置驗(yàn)證

Jezebel-239Pu[10]為239Pu裸球基準(zhǔn)裝置，半徑為6.384 93 cm，基準(zhǔn)實(shí)驗(yàn)的keff為1.000±0.002，裝置核素種類(lèi)及核子密度列于表1。

表1 Jezebel-239Pu裝置核子密度

因NAS系統(tǒng)基于六角形粗網(wǎng)格節(jié)塊法，無(wú)法用于球形裝置計(jì)算，故在Jezebel-239Pu計(jì)算中只使用CITATION程序。計(jì)算中截面數(shù)據(jù)采用33群截面庫(kù)及相應(yīng)的協(xié)方差數(shù)據(jù)，經(jīng)PASC程序[11]計(jì)算共振自屏處理后，得到適合CITATION的輸入格式。

對(duì)于Jezebel-239Pu這種小型裸球基準(zhǔn)裝置，因CITATION基于擴(kuò)散理論，故keff計(jì)算值偏小，利用該套多群截面數(shù)據(jù)計(jì)算得到keff=0.955 2。對(duì)239Pu的輻射俘獲截面σ(n,γ)、σf、σtr進(jìn)行不確定性分析，所得keff分布及正態(tài)擬合結(jié)果示于圖2。由圖2可見(jiàn)，計(jì)算結(jié)果與正態(tài)分布曲線擬合較好。計(jì)算結(jié)果的統(tǒng)計(jì)列于表2(計(jì)算條件為Windows7x64，IVF2011，Intel Core i7-2600，所用程序皆為串行程序)，表2同時(shí)列出了利用S-U方法[12]所得的不確定數(shù)據(jù)作為對(duì)比，表中統(tǒng)計(jì)keff為10 000次計(jì)算結(jié)果的平均值。

圖2 CITATION計(jì)算Jezebel-239Pu裝置keff分布

分析keff計(jì)算不確定度結(jié)果可知，該方法結(jié)果與S-U方法的差別在20%以?xún)?nèi)，差別是因?yàn)橛?jì)算所用截面與協(xié)方差數(shù)據(jù)不同，且S-U方法總不確定度計(jì)算中考慮了其他一些因素引起的不確定度，如各反應(yīng)道之間以及裂變中子數(shù)引入的不確定度，因此所得結(jié)果較本方法的偏大。同時(shí)，在進(jìn)行10 000次計(jì)算條件下，程序所需時(shí)間小于7 min，時(shí)間效率較高。該方法進(jìn)行1個(gè)截面抽樣與進(jìn)行3個(gè)截面抽樣所用時(shí)間一致，因此該方法計(jì)算時(shí)間主要消耗在堆芯計(jì)算程序上，截面抽樣帶來(lái)的額外計(jì)算量很小?？煽吹?，在95%置信區(qū)間下，不確定度范圍很小，所得結(jié)果的精度完全滿(mǎn)足要求。

表2 Jezebel-239Pu裝置keff不確定度

3.2 中國(guó)實(shí)驗(yàn)快堆首爐堆芯驗(yàn)證

中國(guó)實(shí)驗(yàn)快堆首爐堆芯[13]采用富集度64.4%的UO2陶瓷體作為燃料，冷卻劑為液態(tài)金屬鈉，包殼材料為不銹鋼，堆芯裝載示于圖3。

圖3 CEFR初始堆芯裝載布置

計(jì)算中采用的協(xié)方差數(shù)據(jù)為中國(guó)核數(shù)據(jù)中心為CEFR制作的6群截面數(shù)據(jù)，本文對(duì)235U與238U的(n,f)、(n,γ)截面以及56Fe和23Na的(n,γ)截面進(jìn)行抽樣，由于現(xiàn)有數(shù)據(jù)中無(wú)238U的截面協(xié)方差數(shù)據(jù)，故采用與235U相同的相對(duì)協(xié)方差數(shù)據(jù)作為近似。

計(jì)算中截面采用171群NVITAMIN-C庫(kù)，利用PASC系統(tǒng)并群制作適合CEFR的6群截面數(shù)據(jù)。圖4分別示出了NAS和CITATION程序堆芯計(jì)算keff的分布。其中，利用NAS進(jìn)行了10 000次堆芯計(jì)算，利用CITATION進(jìn)行了1 000次堆芯計(jì)算。由圖4可知，兩者的分布均基本符合正態(tài)，但NAS由于統(tǒng)計(jì)次數(shù)多，故符合得更好。

表3列出了CEFR核截面數(shù)據(jù)不確定性造成的keff不確定度結(jié)果，表中標(biāo)準(zhǔn)計(jì)算keff為采用原始數(shù)據(jù)進(jìn)行單次計(jì)算所得結(jié)果。由核數(shù)據(jù)不確定性造成的不確定度為2.38%，作為比照，給出了三維蒙特卡羅微擾方法計(jì)算不確定度為2.27%，俄羅斯設(shè)計(jì)計(jì)算給出的keff不確定度為1.90%[14]，一維S-U方法計(jì)算不確定度為2.65%[15]。計(jì)算結(jié)果與參考值符合較好，其中，文獻(xiàn)[15]及三維蒙特卡羅微擾方法所采用協(xié)方差數(shù)據(jù)與本文一致，可見(jiàn)，該方法所得結(jié)果與參考值中較精確的三維蒙特卡羅微擾方法的結(jié)果的相對(duì)偏差小于10%，好于一維S-U方法的。同時(shí)，NAS與CITATION兩程序計(jì)算結(jié)果很接近，說(shuō)明直接蒙特卡羅方法可用于不同的中子學(xué)計(jì)算程序。

圖4 NAS和CITATION計(jì)算CEFR的keff分布

表3CEFRkeff不確定度計(jì)算結(jié)果

Table3UncertaintyofcalculatedkeffforCEFR

計(jì)算程序統(tǒng)計(jì)keff標(biāo)準(zhǔn)計(jì)算keff相對(duì)不確定度/%95%置信度不確定度范圍/%計(jì)算次數(shù)NAS1.025 81.025 92.38[2.35,2.41]10 000CITATION1.024 51.023 22.35[2.25,2.45]1 000

對(duì)于CITATION程序1 000次計(jì)算，所得95%置信區(qū)間不確定度范圍為[2.25%,2.45%]，對(duì)于NAS程序10 000計(jì)算，所得95%置信區(qū)間不確定度范圍為[2.35%,2.41%]，所得結(jié)果相符較好，且可滿(mǎn)足精度要求，若需進(jìn)一步減小范圍，通過(guò)增加計(jì)算次數(shù)即可達(dá)到。按照蒙特卡羅統(tǒng)計(jì)理論，置信度區(qū)間的寬度與計(jì)算次數(shù)的0.5次方成反比[16]。由于NAS計(jì)算速度較快，故在進(jìn)一步的詳細(xì)計(jì)算中，只采用NAS程序。

表4列出了通過(guò)NAS計(jì)算各主要反應(yīng)道不確定性造成的CEFRkeff不確定度及對(duì)比解，每種情況的計(jì)算次數(shù)為1 000次?？芍?，影響CEFRkeff不確定度的主要截面為235U的(n,γ)截面，其次為235U的(n,f)截面以及238U的(n,γ)截面，堆芯中核子份額較高的Na與Fe的核截面數(shù)據(jù)對(duì)整體keff不確定度的貢獻(xiàn)很小。將計(jì)算所得6個(gè)截面的不確定度利用誤差傳遞公式，平方相加，得出上述截面造成的總不確定度為2.39%，與直接計(jì)算總不確定度2.38%接近。這也驗(yàn)證了該方法既可計(jì)算單個(gè)截面造成的不確定度，也可計(jì)算多個(gè)截面造成的總不確定度。與三維蒙特卡羅微擾方法相比，結(jié)果符合較好，同時(shí)該方法能克服蒙特卡羅微擾方法存在的計(jì)算小的不確定度時(shí)統(tǒng)計(jì)誤差偏大的缺點(diǎn)。

表4 NAS計(jì)算分截面核數(shù)據(jù)不確定性對(duì)CEFR keff不確定度影響

注：1) 采用MCNP微擾方法計(jì)算得到[17]，擾動(dòng)量為10%

4 結(jié)論

核截面數(shù)據(jù)不確定性會(huì)引入核裝置keff計(jì)算值的不確定度。直接蒙特卡羅方法可有效、簡(jiǎn)便地給出所需keff計(jì)算值的不確定度，結(jié)果合理可信。該方法的優(yōu)點(diǎn)在于對(duì)各種裝置keff求解的適應(yīng)性強(qiáng)、可移植性好并可計(jì)算三維問(wèn)題。但與所有蒙特卡羅方法相同，為達(dá)到要求精度，其計(jì)算時(shí)間長(zhǎng)，且存在一定的隨機(jī)性。隨著計(jì)算機(jī)計(jì)算能力以及計(jì)算方法的不斷發(fā)展，該方法的適用范圍會(huì)更加廣泛。

參考文獻(xiàn)：

[1]ALIBERTI G, PALMIOTTI G, SALVATORES M, et al. Nuclear data sensitivity, uncertainty and target accuracy assessment for future nuclear systems[J]. Annals of Nuclear Energy, 2006, 33: 700-733.

[2]許淑艷. 蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用[M]. 北京：原子能出版社，2006.

[3]CHADWICK M B, OBLOZINSKY P, HERMAN M, et al. ENDF/B-Ⅶ.0: Next generation evaluated nuclear data library for nuclear science and technology[J]. Nuclear Data Sheets, 2006, 107(12): 2 931-3 060.

[4]BOX G, MULLER M. A note on the generation of random normal deviates[J]. The Annals of Mathematical Statistics, 1958, 29(2): 610-611.

[5]頁(yè)爾驊，張德平. 概率論與隨機(jī)過(guò)程[M]. 北京：科學(xué)出版社，2005.

[6]茆詩(shī)松，程依明，濮曉龍. 概率論與數(shù)理統(tǒng)計(jì)教程[M]. 北京：高等教育出版社，2004.

[7]李慶揚(yáng)，王能超，易大義. 數(shù)值分析[M]. 5版. 北京：清華大學(xué)出版社，2008.

[8]張棣芳. CITATION程序說(shuō)明書(shū)與使用說(shuō)明[M]. 北京：核工業(yè)核電軟件中心，1995.

[9]李澤華. 三維節(jié)塊法程序在快堆物理計(jì)算中的應(yīng)用[C]∥第八屆反應(yīng)堆數(shù)值與粒子輸運(yùn)學(xué)術(shù)會(huì)議. 深圳：[出版者不詳]，2000.

[10] International handbook of evaluated criticality safety benchmark experiments, NEA/NSC/DOC(95)03[M].US: Nuclear Energy Agency Nuclear Science Committee of the Organization for Economic Cooperation and Development, 2006.

[11] 方邦城，唐忠樑，趙金坤. PASC-1程序系統(tǒng)各模塊功能介紹及其輸入(卡)中文說(shuō)明合訂本，CEFRRLP1.2-CX/PASC-1[R]. 北京：中國(guó)原子能科學(xué)研究院，2006.

[12] 胡澤華，王佳，孫偉力，等. 基準(zhǔn)模型keff對(duì)核數(shù)據(jù)的靈敏度分析及不確定度量化[J]. 原子能科學(xué)技術(shù)，2013，47(增刊):312-317.

HU Zehua, WANG Jia, SUN Weili, et al. Sensitivity and uncertainty analysis of calculatedkeffon benchmark models due to uncertainties of nuclear data[J]. Atomic Energy Science and Technology, 2013, 47(Suppl.): 312-317(in Chinese).

[13] 田和春. 中國(guó)實(shí)驗(yàn)快堆堆芯設(shè)計(jì)說(shuō)明，CEFR01Z19LWS01-SM[R]. 北京：中國(guó)原子能科學(xué)研究院，2002.

[14] 庫(kù)扎夫考夫 Н Г，法拉克欣 М Р，拉玖內(nèi)切娃 А А，等. CEFR堆芯物理特性計(jì)算誤差分析[R]. 俄羅斯：[出版者不詳]，2002.

[15] 剛直. 核截面引起積分參數(shù)keff不確定度的一維分析程序開(kāi)發(fā)[D]. 北京：中國(guó)原子能科學(xué)研究院，2006.

[16] ROBERT C P, CASELLA G. Monte Carlo statistical methods[M]. 2nd ed. New York: Springer, 2004.

[17] 劉萍. 核數(shù)據(jù)不確定性對(duì)ADS系統(tǒng)keff等積分量的影響和WIMS 82群庫(kù)的研制[D]. 北京：中國(guó)原子能科學(xué)研究院，2004.