邱有恒，史 濤

(1.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100094;2.上海核工程研究設(shè)計(jì)院有限公司，上海 200233)

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

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

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

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

事實(shí)上，利用正向蒙特卡羅計(jì)算通量可產(chǎn)生重要性函數(shù)，這方面的優(yōu)秀代表是MCNP程序的權(quán)窗發(fā)生器。本文在MCNP程序基礎(chǔ)上，開發(fā)基于正向蒙特卡羅計(jì)算自動源偏倚抽樣功能，以實(shí)現(xiàn)與權(quán)窗自洽耦合。

1 理論與方法

1.1 偏倚抽樣基本理論

考慮如下數(shù)學(xué)期望的計(jì)算：

(1)

其中：f(x)為分布函數(shù)；g(x)為響應(yīng)函數(shù)，當(dāng)從f(x)抽樣困難或計(jì)算效率很低時(shí)，可考慮如下形式的計(jì)算：

(2)

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

1.2 權(quán)窗與權(quán)窗發(fā)生器

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

權(quán)窗有3個(gè)參數(shù)，分別是權(quán)窗下限(WWL)、權(quán)窗上限(WWU)、存活粒子權(quán)重(WS)。通過權(quán)窗游戲，可將粒子權(quán)重分布控制在合理的范圍內(nèi)。假定某相空間粒子權(quán)重W落在以下3個(gè)區(qū)間，根據(jù)權(quán)窗規(guī)則，粒子權(quán)重將作相應(yīng)調(diào)整。

1)W

2) WWL

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

MCNP程序有權(quán)窗發(fā)生器功能，可自動生成各相空間經(jīng)過優(yōu)化后的權(quán)窗系數(shù)，包含以下兩步。

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

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

(3)

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

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

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

(4)

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

(5)

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

(6)

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

2 數(shù)值算例

2.1 鉛屏蔽模型

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

圖1 鉛屏蔽模型

在相同的計(jì)算機(jī)上，不同模型計(jì)算結(jié)果列于表1，在下文的敘述中，將不采用重要性的方法稱為方法1，將采用幾何重要性的方法稱為方法2，將采用MCNP模式的權(quán)窗方法稱為方法3，將本文方法稱為方法4。采用方法1計(jì)算時(shí)，可同時(shí)產(chǎn)生方法3和方法4的參數(shù)。其中計(jì)算優(yōu)度FOM定義為統(tǒng)計(jì)誤差平方與計(jì)算時(shí)間乘積的倒數(shù)，顯然，F(xiàn)OM越大表明計(jì)算效率越高，表中FOM1表示單次計(jì)算的FOM，F(xiàn)OM2表示累加了產(chǎn)生權(quán)窗參數(shù)所用時(shí)間的FOM。

表1 不同方法的計(jì)算結(jié)果與效率比較

從表1可見，若采用方法1，模擬8×107樣本，統(tǒng)計(jì)誤差較大，達(dá)0.106 3；采用方法2，同樣8×107樣本，統(tǒng)計(jì)誤差很小，僅0.005 7，結(jié)果已收斂，但計(jì)算時(shí)間超長，約21 h，F(xiàn)OM較無重要性的高約30倍。

采用方法3，樣本數(shù)為2×107的情況下，統(tǒng)計(jì)誤差僅0.003 9，低于幾何重要性函數(shù)，計(jì)算結(jié)果與幾何重要性幾乎完全相同，F(xiàn)OM遠(yuǎn)高于幾何重要性。

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

2.2 反應(yīng)堆屏蔽計(jì)算算例

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

圖2 某壓水堆一次屏蔽計(jì)算模型

圖3 中平面處源中子分布

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

圖4 中平面測點(diǎn)的源中子偏倚抽樣概率

圖5 堆芯頂端測點(diǎn)的源中子偏倚抽樣概率

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

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

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

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

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

表2 輻照監(jiān)督管頂端測點(diǎn)通量密度比較

3 結(jié)論

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