蒙特卡羅臨界計算能量偏倚的最佳源偏倚方法

潘清泉，饒俊杰，王 侃

(清華大學 工程物理系，北京 100084)

反應堆臨界問題[1]以求解keff為特征，即求解中子輸運方程的特征值。反應堆理論中，keff被看作相鄰計算代的子代中子數(shù)和父代中子數(shù)的比值，其中不同計算代以裂變反應區(qū)分。蒙特卡羅(MC)臨界計算采用裂變源迭代法[2]，用當前代產(chǎn)生的裂變源作為下一代的初始源，即下一代各裂變源區(qū)抽樣產(chǎn)生的源粒子數(shù)與當前代在該裂變源區(qū)產(chǎn)生的裂變中子數(shù)正相關(guān)，可表示為ni=Npi，N為每個計算代的總中子數(shù)，ni為某計算代第i號裂變源區(qū)所需抽樣產(chǎn)生的源中子數(shù)，pi為第i號裂變源區(qū)產(chǎn)生源中子的概率。如果同時考慮源中子的空間分布和能量分布，則ni,E=Npi,E，ni，E為某計算代第i號裂變源區(qū)抽樣產(chǎn)生的能量為E的源中子數(shù)，pi，E為第i號裂變源區(qū)產(chǎn)生能量為E的源中子的概率。可見，傳統(tǒng)的臨界計算在抽樣每個計算代的源中子時，在空間分布和能量分布上均采用了比例抽樣法。該抽樣法能反映真實的中子輸運物理過程，但計算所得的總體方差在數(shù)學上并未達到最優(yōu)化。同時，物理上的連貫性必定帶來數(shù)學上的相關(guān)性，故MC計算代之間存在較大的相關(guān)性，使得計數(shù)器結(jié)果和keff的估計值存在方差低估計現(xiàn)象[3]。

為實現(xiàn)源中子在空間和能量上的最優(yōu)化分布，消除方差低估計現(xiàn)象，并減少計算方差，計算中主要存在兩個問題：1) 在MC臨界計算的每一計算代中，在某一裂變源區(qū)某一能量區(qū)間中應該抽樣產(chǎn)生多少源中子才能使得總體方差最??；2) 在MC臨界計算的每一計算代中，源中子在空間和能量上怎樣分布才能實現(xiàn)數(shù)學最優(yōu)化。

為此，本文研究MC臨界計算的能量偏倚的最佳源偏倚方法。該方法基于傳統(tǒng)的最佳源偏倚方法，以香農(nóng)熵診斷[4]和組統(tǒng)計[5]方法為基礎(chǔ)，在消除計算代之間相關(guān)性的同時，通過最佳分層抽樣法[6]實現(xiàn)堆用MC程序RMC臨界計算時，每個計算代的源中子在空間和能量上的最佳分布，從而達到全局減方差[7]的計算效果。

1 理論模型

在數(shù)理統(tǒng)計學[8]的隨機抽樣法中，分層抽樣法得到了廣泛關(guān)注。當總體樣本由差異明顯的幾部分組成時，為使樣本更客觀地反映總體情況，通常將總體按照不同的特點分成層次分明的幾部分，然后按照各部分在總體中所占的比例進行抽樣，這種抽樣法稱為分層抽樣法。分層抽樣法是通過改變概率分布，把整層抽樣變?yōu)榉謱映闃?，把整層概率分布變?yōu)榉謱痈怕史植迹瑥亩_到降低方差的計算效果。

傳統(tǒng)的MC方法進行臨界計算時采用了裂變源迭代法，用當前代產(chǎn)生的裂變中子源作為下一代的初始源，即在空間和能量上進行了分層抽樣，并采用了比例分層抽樣法。通過拉格朗日乘子法[9]或Schwarz不等式[9]可推導出最佳分層抽樣法，使計數(shù)器的總體方差最小。

但是，在RMC程序中使用最佳分層抽樣法時需提前知道待求的方差信息，導致方法失效，所以最佳分層抽樣法就像是零方差理論，只能無限接近，并不能真正實現(xiàn)[10]。在保證計算代之間不相關(guān)的前提下，用上一代的方差信息近似當前代的方差信息就可解決這一問題。為實現(xiàn)計算代之間不相關(guān)，需在RMC程序臨界計算過程中消除方差低估計現(xiàn)象。本文從通用性和可靠性的角度出發(fā)，選擇了組統(tǒng)計方法，并根據(jù)香農(nóng)熵診斷確定出了簡單有效的組統(tǒng)計長度。

1.1 最佳分層抽樣法

系統(tǒng)特征值keff及其標準偏差在MC臨界計算的所有活躍代中進行統(tǒng)計，并對所有活躍代的結(jié)果求平均值。keff有3種估計方法，分別為碰撞估計法、吸收估計法和徑跡長度估計法，3種估計值通過統(tǒng)計相關(guān)性進行組合得到無偏的keff及其標準偏差。

設f(x)為區(qū)域D上的概率密度函數(shù)，考慮下列數(shù)學期望的計算：

(1)

其中：R為響應函數(shù)乘以概率密度函數(shù)的積分值；g(x)為響應函數(shù)；x為某一樣本值；E[g]為數(shù)學期望。

(2)

(3)

(4)

這樣把區(qū)域D上的積分轉(zhuǎn)化為m個子區(qū)域的積分之和。

(5)

(6)

(7)

在每一計算代中，由于Ni的計算表達式中涉及待求量σi，因此式(7)很難被直接應用。但使用RMC程序進行臨界計算時，通過非活躍代的計算完成源收斂后，當前代的σi可由上一代的σi近似，故在臨界計算過程中，可采用這種方法實現(xiàn)最佳分層抽樣。即每個計算代在進行裂變源抽樣時可采用最佳源偏倚因子對源中子進行偏倚，由比例抽樣：

ni=piN

(8)

變成最佳分層抽樣：

(9)

相應的最佳源偏倚因子εi為：

(10)

為同時對源中子的空間位置和能量進行基于最佳分層抽樣法的源偏倚，相應的能量偏倚的最佳源偏倚因子變?yōu)椋?/p>

(11)

其中：εi,E′為某計算代第i號裂變源區(qū)能量為E′的最佳源偏倚因子；σi，E′為該計算代第i號裂變源區(qū)能量為E′的計算方差；pj，e為該計算代第i號裂變源區(qū)產(chǎn)生能量為e的源中子的真實物理概率。

從式(11)可知，統(tǒng)計方差大的區(qū)域會被多抽樣，統(tǒng)計方差小的區(qū)域會被少抽樣。因此基于最佳分層抽樣法的源偏倚方法不僅可實現(xiàn)全局減方差的計算效果，還能在一定程度上展平全堆方差的分布。但上述關(guān)于分層抽樣法推導的前提是每個計算代的中子數(shù)足夠大，同時計算代之間相互獨立，故為了實現(xiàn)能量偏倚的最佳源偏倚方法，需采用組統(tǒng)計方法來消除方差低估計的現(xiàn)象。

1.2 組統(tǒng)計方法

(12)

(13)

(14)

(15)

(16)

其中，CR[l]稱為l代滯后協(xié)方差：

CR[l]=E[(xi-E[xi])(xi+l-E[xi+l])]=

cov[xi,xi+l]

(17)

其中，cov[xi,xi+l]為第xi和第xi+l代之間的協(xié)方差。

可見，方差低估計僅取決于活躍代間的間距l(xiāng)，與活躍代位置i無關(guān)。在MC源迭代過程中，每代裂變源分布總是來自于上一代的裂變源分布，即裂變源迭代之間存在正相關(guān)性，因此CR[l]的符號為正，于是得到：

(18)

一般而言，方差低估計問題對keff影響較小，但對通量分布和功率分布等局部計數(shù)器影響較大，這是因為裂變源分布的相關(guān)性對計數(shù)器的影響更為直接。系統(tǒng)占優(yōu)比越高，則源迭代之間的相關(guān)性越強，從而導致方差低估計現(xiàn)象越顯著。

Gelbard等[5]最早提出使用組統(tǒng)計方法來減少MC方差低估計，MC21程序采用了該方法。組統(tǒng)計方法的基本思路是將多個活躍代歸并為1組，然后以組為單位統(tǒng)計keff或計數(shù)器的均值和方差，其具體實現(xiàn)過程如下。

首先假設MC計算共有N個活躍代，按照每組M代進行歸并，則共得到N′=N/M組計數(shù)。每組計數(shù)可表示為：

(19)

其中：xi為第i個活躍代計數(shù)；yj為第j組計數(shù)。

然后，以組為單位統(tǒng)計物理量的平均值和方差：

(20)

(21)

顯然，裂變源分布在組之間的相關(guān)性弱于其在代之間的相關(guān)性，因此組統(tǒng)計方法可減輕裂變源相關(guān)性所引起的方差低估計問題。當每組粒子代數(shù)M足夠大時，方差低估計偏差可忽略不計。

在組統(tǒng)計方法中，關(guān)鍵是確定每組的粒子代數(shù)M，即組長度。一般而言，源迭代過程中裂變源分布的相關(guān)性越強，則意味著源分布需經(jīng)歷更多的代數(shù)才能消除之前的源分布的影響，與之相對應地，所需的組長度M越大。通過計算占優(yōu)比，可定性判斷源分布相關(guān)性的強弱，對組長度M的選擇有一定指導意義。在RMC程序中實現(xiàn)的是另一種很簡便的方法：

M≈Ninact

(22)

其中，Ninact為源收斂所需的非活躍代數(shù)。

MC源收斂速度和方差低估計程度均是由源迭代裂變源分布的相關(guān)性決定的。裂變源的相關(guān)性越強，則源收斂速度越緩慢，同時方差低估計現(xiàn)象也越明顯。MC非活躍代的源收斂過程意味著初始源分布與真實源分布的偏差通過M≈Ninact代的誤差傳遞后被基本消除。同樣地，在MC活躍代中，某一裂變源分布的影響在經(jīng)歷M≈Ninact代后可認為忽略不計。直接使用非活躍代數(shù)確定組統(tǒng)計方法中的組長度M，是符合MC計算的實際情況的。

1.3 香農(nóng)熵診斷

在反應堆物理歷史上很長一段時間內(nèi)占優(yōu)比作為計算收斂速度的重要判斷標準，但后期研究表明收斂的判斷除keff的收斂判斷外，還需包括裂變源分布的收斂判斷，且實際上裂變源分布的收斂要慢于keff，即當裂變源收斂時keff一定收斂，但keff收斂還不能判斷源迭代是否己達到收斂。對于MC臨界計算，應評估keff和裂變源的收斂性。2005年，Ueki等[4]將信息學領(lǐng)域中用于衡量信息量的香農(nóng)熵概念和方法引入到MC粒子輸運方法中，作為臨界計算收斂診斷的方法，當熵收斂后，源分布也就收斂。香農(nóng)熵目前已作為一種具有顯著優(yōu)勢的源收斂診斷方法，在RMC等程序中得到應用。

香農(nóng)熵源收斂診斷方法的核心思想是將每個迭代周期中裂變源中子分布信息用香農(nóng)熵進行宏觀度量，通過觀察臨界迭代計算中香農(nóng)熵隨代數(shù)變化情況來方便地觀察裂變源分布的收斂趨勢，當源分布趨于平穩(wěn)時，香農(nóng)熵收斂于一穩(wěn)態(tài)值。裂變源分布的香農(nóng)熵定義為：

(23)

其中：H(si)為香農(nóng)熵；si為第i代的裂變源分布；B為空間網(wǎng)格劃分的編號。

2 計算流程

能量偏倚的最佳源偏倚方法建立在香農(nóng)熵診斷和組統(tǒng)計方法的基礎(chǔ)上。對于某一特定算例，首先需通過香農(nóng)熵診斷確定出組統(tǒng)計的長度M，然后在MC臨界計算中設置組統(tǒng)計參數(shù)。在裂變源迭代法計算過程中，每完成1次組統(tǒng)計就會得到各裂變源區(qū)在各能量區(qū)間內(nèi)的方差信息σi,E。然后，所需的當前代方差信息由上一代的方差信息近似得到，根據(jù)式(11)計算出各裂變源區(qū)在各能量區(qū)間內(nèi)的最佳源偏倚因子。在保證每個計算代中子總權(quán)重守恒的基礎(chǔ)上，對抽樣得到的裂變源中子在空間和能量上進行偏倚。假設某個裂變源區(qū)在某能量區(qū)間內(nèi)的最佳源偏倚因子為εi,E，則可通過簡單的輪盤賭與分裂實現(xiàn)最佳源偏倚。

(24)

其中，ωbias和ωbank分別為源偏倚后的粒子權(quán)重和存庫的粒子權(quán)重。

以上的處理方式并不會改變每個計算代裂變中子的總權(quán)重，但實現(xiàn)了源中子在空間上的最佳分布。具體的計算流程如下：1) 在RMC程序的輸入文件中設置能量分箱；2) 采用香農(nóng)熵診斷確定組長度M；3) 設置好組統(tǒng)計的參數(shù)，開始臨界計算；4) 用組統(tǒng)計的方差信息計算最佳源偏倚因子；5) 對粒子的空間位置和能量進行偏倚；6) 回到第4步，更新最佳源偏倚因子；7) 迭代計算，直到計算收斂。能量偏倚的最佳源偏倚方法的計算流程如圖1所示。

圖1 能量偏倚的最佳源偏倚方法的計算流程

3 數(shù)值驗證

選擇典型17×17的壓水堆組件[12]進行臨界計算。燃料棒是密度為10.196 g/cm3的UO2，鋯合金的密度為6.550 g/cm3。幾何建模采用17×17的柵元劃分，共有289個網(wǎng)格。為更好地展示能量偏倚的最佳源偏倚方法對方差的展平功能，邊界條件選擇了真空邊界條件。組件的幾何結(jié)構(gòu)如圖2所示。

圖2 典型壓水堆組件的幾何結(jié)構(gòu)

首先對粒子進行能量分箱，本次計算采用C5G7模型[13]的能量分箱，共劃分了7個能量區(qū)間(單位MeV)，從小到大依次為：[0.0，1.3×10-7]，[1.3×10-7，6.3×10-7]，[6.3×10-7，4.1×10-6]，[4.1×10-6，5.56×10-5]，[5.56×10-5，9.2×10-3]，[9.2×10-3，1.36]，[1.36，10]。對該組件進行香農(nóng)熵統(tǒng)計，每代50 000個源中子，共進行了550個非活躍代的香農(nóng)熵統(tǒng)計，計算結(jié)果如圖3所示。由圖3可看出，經(jīng)過十幾次非活躍代的迭代過程，組件的源分布基本達到收斂，故組統(tǒng)計方法選擇了組長度M=20。

然后使用RMC程序?qū)υ摻M件進行能量偏倚的最佳源偏倚方法進行計算，計算條件為50個非活躍代，2 010個活躍代，其中組統(tǒng)計的長度M=20，所以程序進行了98次組統(tǒng)計計算。本文同時也進行了每代統(tǒng)計(M=1)的最佳源偏倚方法計算。作為對比，使用RMC程序進行了比例抽樣法的臨界計算，也即未進行最佳源偏倚方法的計算。

圖3 典型壓水堆組件香農(nóng)熵統(tǒng)計的計算結(jié)果

能量偏倚的最佳源偏倚方法可實現(xiàn)全局減方差的計算效果，為更好地量化描述全局減方差的效果，本文選擇如下參考量[11]。

平均相對偏差：

(25)

平均品質(zhì)因子：

(26)

相對偏差標準差：

(27)

其中：N為計算結(jié)果的個數(shù)；Rei為第i個計算結(jié)果的統(tǒng)計相對偏差；T為計算時間。

表1列出壓水堆組件能量偏倚的最佳源偏倚方法的計算結(jié)果。由表1可得出如下結(jié)論。

表1 壓水堆組件能量偏倚的最佳源偏倚方法的計算結(jié)果

1) 不論是每代統(tǒng)計還是組統(tǒng)計，能量偏倚的最佳源偏倚方法均能保證計算結(jié)果的精度，即能量偏倚的最佳源偏倚方法的計算結(jié)果與未進行能量偏倚的最佳源偏倚的計算結(jié)果相同。

2) 不論是每代統(tǒng)計還是組統(tǒng)計，能量偏倚的最佳源偏倚方法均能展平方差的分布。因為能量偏倚的最佳源偏倚因子使得方差大的裂變源區(qū)多抽樣，方差小的裂變源區(qū)少抽樣。

3) 基于每代統(tǒng)計的能量偏倚的最佳源偏倚方法與比例抽樣法的平均相對偏差、平均品質(zhì)因子和相對偏差標準差各有大小，未見明顯區(qū)別。因為每代統(tǒng)計的計算結(jié)果存在方差低估計現(xiàn)象，代與代之間具有較大相關(guān)性。

4) 基于組統(tǒng)計的能量偏倚的最佳源偏倚方法，平均相對偏差、平均品質(zhì)因子和相對偏差標準差明顯優(yōu)于比例抽樣法，可見消除了代與代之間的相關(guān)性后，能量偏倚的最佳源偏倚方法可有效實現(xiàn)全局減方差的效果。

4 結(jié)論

使用MC程序進行反應堆臨界計算時，裂變源迭代法會帶來方差低估計的現(xiàn)象。通過分析源中子的抽樣過程，以及對分層抽樣法進行數(shù)學推導，本文建立起了能量偏倚的最佳源偏倚方法。該方法基于香農(nóng)熵診斷和組統(tǒng)計方法，對源中子的空間位置和能量進行偏倚，展平了全堆方差分布，并實現(xiàn)了全局減方差的計算效果。對一標準壓水堆組件進行了測試，計算結(jié)果證明了能量偏倚的最佳源偏倚方法的正確性和有效性。