基于MESH技術(shù)的伴隨權(quán)重窗自動(dòng)生成方法及其應(yīng)用

邱有恒，鄧 力，李百文，趙英奎

(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所 計(jì)算物理實(shí)驗(yàn)室，北京 100094)

對(duì)于大型復(fù)雜系統(tǒng)或深穿透問(wèn)題，直接蒙特卡羅方法的計(jì)算效率很低。通過(guò)設(shè)置合理的重要性函數(shù)，用于指導(dǎo)粒子的賭與分裂，使重要的粒子多抽樣，不重要的粒子少抽樣，可大幅提升計(jì)算效率[1-2]。

重要性函數(shù)反映從相空間(包括幾何空間、能量、時(shí)間等)發(fā)射的粒子的相對(duì)重要性，其實(shí)質(zhì)是偏倚函數(shù)。顯然，理想的重要性函數(shù)值應(yīng)是正比于該相空間內(nèi)的粒子對(duì)探測(cè)器的貢獻(xiàn)。但在獲得精確的數(shù)值模擬結(jié)果前，各相空間的重要性函數(shù)值未知，可通過(guò)多次迭代逐步逼近獲得理想的重要性函數(shù)和精確數(shù)值模擬結(jié)果[3-4]。

幾何重要性函數(shù)和權(quán)重窗技巧[5]是兩種常用的重要性函數(shù)。與幾何重要性函數(shù)相比，權(quán)重窗技巧的本領(lǐng)更強(qiáng)，它能反映不同幾何空間、能量間隔、時(shí)間間隔的粒子的重要性，既能作用于柵元界面，也能作用于碰撞點(diǎn)。且MCNP程序有權(quán)重窗發(fā)生器功能，可多次迭代自動(dòng)產(chǎn)生優(yōu)化的權(quán)重窗參數(shù)。因此，在解決大型復(fù)雜系統(tǒng)或深穿透問(wèn)題中，權(quán)重窗技術(shù)較幾何重要性函數(shù)更靈活、方便，效果更好。

MCNP程序有兩種權(quán)重窗發(fā)生器，一種是正算權(quán)重窗發(fā)生器，另一種是多群伴隨權(quán)重窗發(fā)生器。目前，MCNP程序只能提供基于柵元的伴隨權(quán)重窗發(fā)生器。對(duì)一些含重復(fù)結(jié)構(gòu)或非對(duì)稱結(jié)構(gòu)柵元的模型，由于MCNP程序不能自動(dòng)計(jì)算他們的體積，伴隨權(quán)重窗發(fā)生器將失效。

本文借鑒MCNP程序正算MESH權(quán)重窗發(fā)生功能[6]，對(duì)多群伴隨權(quán)重窗發(fā)生器進(jìn)行基于MESH技術(shù)的伴隨權(quán)重窗自動(dòng)產(chǎn)生功能，并應(yīng)用于光子散射基準(zhǔn)模型模擬。

1 權(quán)重窗技術(shù)

權(quán)重窗是諸多降低方差技巧中的一種[4]。權(quán)重窗有權(quán)重窗下限WL、權(quán)重窗上限WU、存活粒子權(quán)重WS3個(gè)參數(shù)。若粒子權(quán)重低于WL，將對(duì)粒子實(shí)施輪盤(pán)賭游戲；若粒子權(quán)重高于上限WU，將對(duì)粒子實(shí)施分裂游戲；其他情況不做處理。權(quán)重窗對(duì)不重要的粒子少花計(jì)算資源，節(jié)省計(jì)算時(shí)間；對(duì)重要粒子多花計(jì)算資源，從而降低計(jì)算方差，并最終提高計(jì)算效率。

2 伴隨權(quán)重窗自動(dòng)產(chǎn)生方法

2.1 伴隨通量作為重要性函數(shù)的理論基礎(chǔ)

蒙特卡羅方法求解粒子輸運(yùn)方程通常是解發(fā)射密度型積分方程[2]：

(1)

其中：ω(s)為發(fā)射密度；K(s′→s)為轉(zhuǎn)移函數(shù)；Q(s)為源。各種響應(yīng)量I(反應(yīng)率、劑量率等)均可通過(guò)響應(yīng)函數(shù)g(s)與發(fā)射密度的卷積實(shí)現(xiàn)，有：

(2)

與之對(duì)應(yīng)的伴隨輸運(yùn)方程為：

(3)

K+(s′→s)=K(s→s′)

(4)

將式(1)乘以ω+(s),式(3)乘以ω(s)，相減后積分，可得：

(5)

可見(jiàn)，伴隨發(fā)射密度ω+(s)的意義是1個(gè)源粒子對(duì)I的貢獻(xiàn)，也可稱為價(jià)值函數(shù)。

伴隨通量與伴隨發(fā)射密度的關(guān)系為：

(6)

其中：T(r′→r|Ω,E)為輸運(yùn)核；Σt(r,E)為總截面。因此，伴隨通量是天然的重要性函數(shù)。

2.2 從伴隨通量產(chǎn)生權(quán)重窗的方法

由于不同柵元、不同能群的伴隨通量可能有量級(jí)的差別，若直接使用伴隨通量指導(dǎo)粒子輸運(yùn)可能造成粒子徑跡或權(quán)重的大幅波動(dòng)。因此，需對(duì)伴隨通量進(jìn)行處理以避免大幅波動(dòng)。首先，對(duì)伴隨通量進(jìn)行一系列平滑處理，使相鄰相空間伴隨通量之比不過(guò)大。然后，根據(jù)相應(yīng)的參數(shù)產(chǎn)生權(quán)重窗下限系數(shù)，基本思想是權(quán)重窗下限反比于伴隨通量。假定參考柵元中最大的伴隨群通量為e，用戶指定的參考柵元權(quán)重窗下限參數(shù)為fnw(默認(rèn)為1)，平滑指數(shù)為pm(0～1)，按式(7)產(chǎn)生柵元i、能群g對(duì)應(yīng)的權(quán)重窗下限參數(shù)。

(7)

參考柵元一般為正算的源柵元，若fnw取默認(rèn)值1，源柵元的權(quán)重窗下限均≥1，即源粒子均將進(jìn)行輪盤(pán)賭，存活下來(lái)的粒子權(quán)重上升，可游動(dòng)到更遠(yuǎn)的地方。

2.3 基于MESH技術(shù)的伴隨權(quán)重窗實(shí)現(xiàn)方法

利用MCNP程序的FMESH功能可得到基于虛擬網(wǎng)格技術(shù)的伴隨通量分布，由此可避開(kāi)重復(fù)結(jié)構(gòu)或其他非對(duì)稱結(jié)構(gòu)不能計(jì)算體積的難點(diǎn)。MESH網(wǎng)格是圓柱體或長(zhǎng)方體，每個(gè)小的虛擬網(wǎng)格為規(guī)則的幾何體，便于計(jì)算體積。對(duì)每個(gè)小的虛擬網(wǎng)格，采用徑跡長(zhǎng)度估計(jì)方法計(jì)算體通量。

在模塊fmesh_mod中增加一子程序amesh，在amesh中對(duì)各相空間伴隨通量進(jìn)行平滑等處理，產(chǎn)生最終的權(quán)重窗下限，并按正算mesh-based權(quán)重窗發(fā)生器的格式輸出在wwout文件中。程序經(jīng)修改后，在輸入卡片inp中只要同時(shí)有正算的mesh權(quán)重窗發(fā)生器參數(shù)和fmesh計(jì)數(shù)，并采用多群伴隨輸運(yùn)模式，即可產(chǎn)生mesh-based伴隨權(quán)重窗。相關(guān)的子程序調(diào)用情況大致如下：

1) 子程序avrwgi輸出mesh權(quán)重窗定義參數(shù)，如能群、網(wǎng)格劃分等，由于輸入卡片中有正算mesh權(quán)重窗產(chǎn)生指令，這一步會(huì)正常進(jìn)行；

2) 在子程序summary中，先調(diào)用子程序fmesh_print輸出fmesh計(jì)數(shù)結(jié)果，然后調(diào)用權(quán)重窗輸出子程序avrout；

3) 在avrout中，若滿足既有mesh權(quán)重窗，又是伴隨輸運(yùn)，則調(diào)用amesh產(chǎn)生伴隨mesh權(quán)重窗系數(shù)；

4) 在amesh中，利用fmesh的多維通量記錄數(shù)組完成伴隨通量平滑并產(chǎn)生最終權(quán)重窗系數(shù)。

3 數(shù)值算例

3.1 重復(fù)結(jié)構(gòu)模型驗(yàn)證

圖1為MCNP程序中自帶的一重復(fù)結(jié)構(gòu)例題，大致結(jié)構(gòu)為半徑15 cm的球，上、下半球分別被兩種格子模型填充，1號(hào)柵元內(nèi)各向同性發(fā)射均勻球體積源，測(cè)點(diǎn)在上半球。

圖1 重復(fù)結(jié)構(gòu)模型

若采用幾何重要性函數(shù)(imp卡片)，則下半球中所有的8號(hào)小球區(qū)具有相同的重要性，所有的9號(hào)柵元(空白區(qū))具有相同重要性，上半球中左側(cè)所有7號(hào)格具有相同重要性，右側(cè)所有7號(hào)格具有相同重要性，這顯然是粗糙甚至不合理的重要性函數(shù)設(shè)置。為便于比較，將各柵元幾何重要性全設(shè)置為1，即等價(jià)于無(wú)重要性函數(shù)。

若采用MCNP程序基于柵元的多群伴隨輸運(yùn)權(quán)重窗功能(cell-based)，則由于不能計(jì)算體積無(wú)法開(kāi)展計(jì)算。由注量率和計(jì)算效率的比較(表1)可知，采用mesh-based伴隨權(quán)重窗后，計(jì)算結(jié)果幾乎不變，統(tǒng)計(jì)誤差更小。為便于比較，常用優(yōu)度FOM[4]反映計(jì)算效率，F(xiàn)OM越大，計(jì)算效率越高。采用伴隨權(quán)重窗后的FOM約為無(wú)重要性函數(shù)的8.77倍。

圖2為其中1個(gè)能群(0.82～1.35 MeV)的伴隨權(quán)重窗下限分布云圖。權(quán)重窗下限系數(shù)反比于重要性，從圖2可見(jiàn)，測(cè)點(diǎn)附近區(qū)域的重要性最高(權(quán)重窗下限最低，分裂概率最大)，距測(cè)點(diǎn)越遠(yuǎn)重要性越低，4個(gè)角上的深色區(qū)域?yàn)槌瞿Ｐ偷牟糠?，重要性?，權(quán)重窗下限也為0。

表1 注量率和優(yōu)度的比較

圖2 能群z=0平面伴隨權(quán)重窗下限分布云圖

3.2 光子空氣散射基準(zhǔn)問(wèn)題

為檢驗(yàn)數(shù)值模擬程序和參數(shù)精度，美國(guó)堪薩斯州立大學(xué)于1977年開(kāi)展了光子在空氣中散射的基準(zhǔn)實(shí)驗(yàn)[7]，實(shí)驗(yàn)采用60Co點(diǎn)源，位于水泥柱殼的對(duì)稱軸上高出地面198.12 cm處，從源出發(fā)的光子被水泥柱殼校準(zhǔn)成150.5°方向上的一圓錐束，水泥柱殼足夠厚以至于穿出柱殼的光子可忽略不計(jì)(圖3)。探測(cè)器分別放在水泥柱殼一側(cè)的50、100、200、300、400、500、600、700 m遠(yuǎn)離地面上方1 m處。文獻(xiàn)[7-8]分別采用不同版本MCNP程序模擬了該問(wèn)題，文獻(xiàn)[7]模擬結(jié)果與實(shí)驗(yàn)結(jié)果符合較好，但未公布FOM。文獻(xiàn)[8]采用并行化的MCNP4C版本，樣本數(shù)為1億，模擬結(jié)果與實(shí)驗(yàn)結(jié)果相差較遠(yuǎn)，具體原因未知。為提高計(jì)算效率，本文考慮人為設(shè)定幾何重要性，最大達(dá)400(圖4)。

圖3 水泥柱殼剖面圖

圖4 幾何重要性設(shè)置

本文取最難計(jì)算點(diǎn)(700 m,0 m,1 m)作為對(duì)比計(jì)算模型，首先以此點(diǎn)作為伴隨輸運(yùn)的源點(diǎn)，并采用角度偏倚抽樣方案，獲得MESH格式下各相空間伴隨通量，據(jù)此生成權(quán)重窗下限系數(shù)，圖5示出了其中較重要能群(0.8～1 MeV)在z=1 m這層幾何上權(quán)重窗下限系數(shù)分布，其中，源點(diǎn)是正算的源點(diǎn)和伴隨計(jì)算的參考點(diǎn)，探測(cè)器是伴隨的源點(diǎn)和正算的測(cè)點(diǎn)?？梢?jiàn)，距離探測(cè)器越近的空間，權(quán)重窗下限系數(shù)越低，即這些區(qū)域更重要。

本文采用MCNP5程序在計(jì)算機(jī)上分別開(kāi)展無(wú)重要性的直接模擬和采用伴隨權(quán)重窗技巧模擬。采用直接模擬，即便樣本數(shù)達(dá)107，MCNP程序規(guī)定的10項(xiàng)統(tǒng)計(jì)指標(biāo)有error、vov、slope 3項(xiàng)未達(dá)要求，計(jì)算結(jié)果和統(tǒng)計(jì)誤差與文獻(xiàn)[8]相當(dāng)(表2)，但FOM大很多，這可能與程序版本有關(guān)。

采用伴隨權(quán)重窗技巧模擬時(shí)，樣本數(shù)僅5×106，10項(xiàng)統(tǒng)計(jì)指標(biāo)全部通過(guò)，說(shuō)明收斂正常，結(jié)果可信度高，伴隨權(quán)重窗結(jié)果與實(shí)驗(yàn)吻合很好(已將文獻(xiàn)[7]中照射率單位轉(zhuǎn)換為國(guó)際單位，從通量到照射率的轉(zhuǎn)換公式可參考文獻(xiàn)[7])，相對(duì)偏差僅1.2%。采用伴隨權(quán)重窗的FOM約為直接模擬的8倍，計(jì)算精度和效率均高于文獻(xiàn)[8]。

圖5 能群z=1 m平面伴隨權(quán)重窗下限分布云圖

表2 本文與文獻(xiàn)的照射率和優(yōu)度比較

4 結(jié)論

本文在MCNP程序上實(shí)現(xiàn)了基于MESH技術(shù)的伴隨權(quán)重窗自動(dòng)產(chǎn)生功能，彌補(bǔ)了MCNP程序多群伴隨權(quán)重窗發(fā)生器不能用于含重復(fù)結(jié)構(gòu)等復(fù)雜幾何模型的缺陷。將基于MESH技術(shù)的伴隨權(quán)重窗技巧應(yīng)用于光子散射基準(zhǔn)實(shí)驗(yàn)?zāi)M，計(jì)算結(jié)果與實(shí)驗(yàn)值吻合很好，相對(duì)偏差僅1.2%，計(jì)算效率約為直接模擬的8倍。

參考文獻(xiàn)：

[1] 謝仲生，鄧力. 中子輸運(yùn)理論數(shù)值計(jì)算方法[M]. 西安：西北工業(yè)大學(xué)出版社，2005.

[2] 杜書(shū)華，等. 輸運(yùn)問(wèn)題的計(jì)算機(jī)模擬[M]. 長(zhǎng)沙：湖南科學(xué)技術(shù)出版社，1988.

[3] MENDIUS P W. Group IS-11, MCNP: Multigroup/adjoint capabilities, LA-12704[R]. USA: Los Alamos National Laboratory, 1994.

[4] BRIESMEISTER J F. MCNP: A general Monte Carlo code forN-particle transport, LA-12625-M[R]. USA: Los Alamos National Laboratory, 2000.

[5] HENDRICKS J S, CULBERTSON C N. An assessment of MCNP weight windows, LA-UR-00-55[R]. USA: Los Alamos National Laboratory, 2000.

[6] LIU L R, GARDNER P. A geometry-indepen-dent fine-mesh-based Monte Carlo importance generator[J]. Nucl Sci Eng, 1997, 125(2): 188-195.

[7] OLSHER R H, HSU H H, HARVEY W F. Benchmarking the MCNP Monte Carlo code with a photon skyshine experiment[J]. Nucl Sci Eng, 1993, 114(3): 219-227.

[8] 潘流俊，鄧力. 光子散射基準(zhǔn)問(wèn)題的Monte Carlo模擬[J]. 清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2007， 47(增刊)：955-959.

PAN Liujun, DENG Li. Monte Carlo simulations of the photon skyshine benchmark problem[J]. J Tsinghua Univ: Sci & Tech, 2007, 47(Suppl.): 955-959(in Chinese).