基于正向蒙特卡羅計算的自動源偏倚方法在屏蔽計算中的應用

邱有恒，史 濤

(1.北京應用物理與計算數學研究所，北京 100094;2.上海核工程研究設計院有限公司，上海 200233)

屏蔽計算問題是蒙特卡羅方法的難點問題。以反應堆屏蔽計算為例，僅有少量處于堆芯外圍組件的高能中子能到達屏蔽層外，若采用無偏抽樣，大量的計算時間用于模擬無用的源粒子，計算效率很低。

偏倚抽樣是提高蒙特卡羅模擬計算效率最核心的思想，包括3類偏倚方案[1]：源偏倚、輸運偏倚(如賭與分裂、指數變換)、碰撞偏倚(如離散角度偏倚、強迫碰撞、隱俘獲等)。

在學術界和核工程領域，專家和學者發展了大量的降低方差技巧，包括指數變換法、分層抽樣法、輪盤賭和分裂方法、統計估計法、權重窗等[2-4]，但這些方法主要集中在輸運和碰撞偏倚。MCNP程序[5]是通用的粒子輸運蒙特卡羅程序，配備了多種降低方差技巧，但同樣缺乏自動源偏倚功能。而對于反應堆屏蔽這樣的問題，僅采用輸運偏倚和碰撞偏倚通常是不夠的，必須要配合源偏倚。雖然用戶也可利用個人經驗進行源偏倚，但缺乏科學量化依據且難以處理復雜源。

對于自動源偏倚抽樣比較有代表性的方法是CADIS方法[6]，該方法是利用確定論方法開展伴隨輸運，利用伴隨通量實現源偏倚抽樣。JMCT軟件也配備了基于伴隨通量的自動源偏倚抽樣技巧[7]，該軟件是用蒙特卡羅方法進行伴隨輸運，采用虛擬MESH方法存放偏倚參數，適用于各種復雜源描述。該軟件在反應堆屏蔽計算中取得了很好的效果。

事實上，利用正向蒙特卡羅計算通量可產生重要性函數，這方面的優秀代表是MCNP程序的權窗發生器。本文在MCNP程序基礎上，開發基于正向蒙特卡羅計算自動源偏倚抽樣功能，以實現與權窗自洽耦合。

1 理論與方法

1.1 偏倚抽樣基本理論

考慮如下數學期望的計算：

(1)

其中：f(x)為分布函數；g(x)為響應函數，當從f(x)抽樣困難或計算效率很低時，可考慮如下形式的計算：

(2)

其中：f′(x)為偏倚分布函數；f(x)/f′(x)為糾偏因子。正確使用偏倚抽樣能保證數學期望不變，計算效率提高。

1.2 權窗與權窗發生器

權窗是諸多降低方差技巧中的一種。權窗可替換幾何重要性和能量分裂來指定空間和能量重要性函數，也可提供時間關聯的重要性函數，較幾何重要性函數的信息量更大。

權窗有3個參數，分別是權窗下限(WWL)、權窗上限(WWU)、存活粒子權重(WS)。通過權窗游戲，可將粒子權重分布控制在合理的范圍內。假定某相空間粒子權重W落在以下3個區間，根據權窗規則，粒子權重將作相應調整。

1)W

2) WWL

3)W>WWU，對粒子實施分裂游戲。首先分裂成[W/WS]個子粒子([ ]表示取整)，再以概率W/WS-[W/WS]產生另一個子粒子，它們的權重均為WS。

MCNP程序有權窗發生器功能，可自動生成各相空間經過優化后的權窗系數，包含以下兩步。

1) 統計各相空間重要性。對于柵元i的第g群粒子，其重要性等于由該相空間發出的所有徑跡對探測器的貢獻除以進入該相空間的所有粒子數，用IMP(r,E)表示。

2) 由重要性轉換為權窗下限。選定1個參考柵元(通常為源所在柵元)，確定參考柵元最大重要性的相空間的權窗下限wg，wg由用戶指定，默認值為源粒子平均權重的1/2。其他所有相空間權窗下限反比于其自身的重要性。

(3)

式中：max(IMP(rref,E))為參考柵元不同能群源粒子對探測器貢獻的最大值；wg默認為0.5，因此，對于源區，權窗下限最小值為0.5。

1.3 基于正向蒙特卡羅計算重要性的源偏倚抽樣方法

本文采用與JMCT和CADIS方法同樣的思想產生源偏倚函數，即某相空間偏倚概率正比于重要性函數與無偏概率乘積，只是此處用正向蒙特卡羅計算過程產生的重要性函數替換伴隨通量。

(4)

其中：S(r,E)為無偏倚的抽樣概率；Sb(r,E)為偏倚抽樣概率。源粒子權重修正：

(5)

修正后的源粒子權重可能很大或很小，必須采用與之匹配的權窗系數，否則將發生大量的賭與分裂，從而失去源偏倚的意義。權窗下限為：

(6)

MCNP程序的Cu默認值為5，即源區權窗下限為糾偏后權重的1/3，確保偏倚抽樣后糾偏后的源粒子不會發生賭；由于MCNP程序默認的權窗上限是下限的5倍，因此也不會發生分裂。對于非源區，權窗下限直接反比于重要性。

2 數值算例

2.1 鉛屏蔽模型

如圖1所示，模型為150 cm高的圓柱，圓柱內為鉛，源點位于圓柱底部半徑25 cm圓盤內，源中子能譜為裂變譜，測點坐標(0，0，155 cm)。

圖1 鉛屏蔽模型

在相同的計算機上，不同模型計算結果列于表1，在下文的敘述中，將不采用重要性的方法稱為方法1，將采用幾何重要性的方法稱為方法2，將采用MCNP模式的權窗方法稱為方法3，將本文方法稱為方法4。采用方法1計算時，可同時產生方法3和方法4的參數。其中計算優度FOM定義為統計誤差平方與計算時間乘積的倒數，顯然，FOM越大表明計算效率越高，表中FOM1表示單次計算的FOM，FOM2表示累加了產生權窗參數所用時間的FOM。

表1 不同方法的計算結果與效率比較

從表1可見，若采用方法1，模擬8×107樣本，統計誤差較大，達0.106 3；采用方法2，同樣8×107樣本，統計誤差很小，僅0.005 7，結果已收斂，但計算時間超長，約21 h，FOM較無重要性的高約30倍。

采用方法3，樣本數為2×107的情況下，統計誤差僅0.003 9，低于幾何重要性函數，計算結果與幾何重要性幾乎完全相同，FOM遠高于幾何重要性。

采用方法4，計算結果與幾何重要性結果幾乎完全一致，計算統計誤差較方法3的更小，但FOM更大，尤其FOM2較方法3約提高2倍。

2.2 反應堆屏蔽計算算例

圖2為某商用壓水堆一次屏蔽計算模型。實驗測點位于30°方向上第6根輻照監督管中平面和頂端。中平面測點的數值模擬相對容易，與測試結果也符合較好，但頂端測點的數值模擬收斂很慢，采用幾何重要性方法模擬2×109樣本(MCNP上限)所得結果與測試值相差較大。為此，本文采用基于正向蒙特卡羅計算重要性函數的自動源偏倚與權窗耦合技巧，計算結果與測試值符合較好，計算精度和效率遠高于幾何重要性方法。

圖2 某壓水堆一次屏蔽計算模型

圖3 中平面處源中子分布

圖3示出了堆芯中平面處pin_by_pin源中子分布。圖4示出了中平面測點的源中子偏倚抽樣概率，從圖中可見，僅少量外圍組件有抽樣，絕大部分組件抽樣概率為0。圖5為頂端測點的源中子抽樣概率分布，與圖4現象類似。

圖4 中平面測點的源中子偏倚抽樣概率

圖5 堆芯頂端測點的源中子偏倚抽樣概率

圖6 堆芯源中子抽樣能譜比較

圖6示出了源中子能量偏倚概率，對于中平面測點，0.4 MeV以下源中子抽樣概率為0；對于頂部測點，1 MeV以下源中子抽樣概率為0。通過能量偏倚抽樣，使更多高能源中子參與輸運，提高計算效率。

圖7示出了堆芯軸向抽樣概率比較。無偏抽樣時，堆芯中部抽樣概率高，頂端和低端抽樣概率低；針對中平面測點的偏倚抽樣主要集中在z為-60～40 cm之間，其余部分抽樣概率幾乎為0；針對頂端的偏倚抽樣集中在z大于80 cm的部分。

圖7 堆芯z軸抽樣概率比較

表2列出了頂端測點通量密度結果比較(中平面測點難度較小，未列出)。若只采用幾何重要性，即便樣本數達到2×109，各能群中子通量密度均明顯小于測試值，其原因是源中子高能段抽樣不夠；采用源偏倚與權窗耦合技巧后，計算效率大幅提高，只需2×107樣本，在單機上耗時約0.5 h，各能群中子通量密度統計誤差均在2%以內，收斂狀況良好，中子通量密度計算值與實驗值相對誤差均小于20%，達到工程設計的要求，計算優度FOM為93，較只采用幾何重要性的計算優度0.68提高2個量級。

表2 輻照監督管頂端測點通量密度比較

3 結論

本文利用MCNP程序的權窗發生器產生的重要性函數，生成了自動源偏倚以及與之耦合的權窗系數，在屏蔽計算中取得了很好的效果。本文方法中偏倚參數既可采用柵元模式也可采用虛擬MESH模式，與MCNP的權窗功能完全兼容，使用方便。