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

潘清泉，饒俊杰，王 侃

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

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

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

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

1 理論模型

在數理統計學[8]的隨機抽樣法中，分層抽樣法得到了廣泛關注。當總體樣本由差異明顯的幾部分組成時，為使樣本更客觀地反映總體情況，通常將總體按照不同的特點分成層次分明的幾部分，然后按照各部分在總體中所占的比例進行抽樣，這種抽樣法稱為分層抽樣法。分層抽樣法是通過改變概率分布，把整層抽樣變為分層抽樣，把整層概率分布變為分層概率分布，從而達到降低方差的計算效果。

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

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

1.1 最佳分層抽樣法

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

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

(1)

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

(2)

(3)

(4)

這樣把區域D上的積分轉化為m個子區域的積分之和。

(5)

(6)

(7)

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

ni=piN

(8)

變成最佳分層抽樣：

(9)

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

(10)

為同時對源中子的空間位置和能量進行基于最佳分層抽樣法的源偏倚，相應的能量偏倚的最佳源偏倚因子變為：

(11)

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

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

1.2 組統計方法

(12)

(13)

(14)

(15)

(16)

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

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

cov[xi,xi+l]

(17)

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

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

(18)

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

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

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

(19)

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

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

(20)

(21)

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

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

M≈Ninact

(22)

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

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

1.3 香農熵診斷

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

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

(23)

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

2 計算流程

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

(24)

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

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

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

3 數值驗證

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

圖2 典型壓水堆組件的幾何結構

首先對粒子進行能量分箱，本次計算采用C5G7模型[13]的能量分箱，共劃分了7個能量區間(單位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]。對該組件進行香農熵統計，每代50 000個源中子，共進行了550個非活躍代的香農熵統計，計算結果如圖3所示。由圖3可看出，經過十幾次非活躍代的迭代過程，組件的源分布基本達到收斂，故組統計方法選擇了組長度M=20。

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

圖3 典型壓水堆組件香農熵統計的計算結果

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

平均相對偏差：

(25)

平均品質因子：

(26)

相對偏差標準差：

(27)

其中：N為計算結果的個數；Rei為第i個計算結果的統計相對偏差；T為計算時間。

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

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

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

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

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

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

4 結論

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