韓靜茹，陳義學(xué)，石生春，袁龍軍，陸道綱

（1．華北電力大學(xué)核科學(xué)與工程學(xué)院，北京 102206；2．環(huán)境保護部核與輻射安全中心，北京100082）

輻射屏蔽系統(tǒng)設(shè)計是核裝置工程設(shè)計的核心內(nèi)容之一，其設(shè)計優(yōu)劣直接影響工程造價及工作人員與周圍環(huán)境的輻射安全。目前常用的輻射屏蔽計算方法分為兩種：確定性方法（如SN）和蒙特卡羅方法（MC）。其中SN方法特別適合求解幾何相對簡單的“深穿透”問題，并可提供整個模型的分布解。然而，其在處理復(fù)雜幾何模擬上存在限制，同時計算時間和所需存儲單元隨維度增加而增加。MC方法的優(yōu)勢在于能夠精確描述復(fù)雜幾何結(jié)構(gòu)，且其計算量不受維度的影響。但由于其抽樣本質(zhì)，對于大型系統(tǒng)深穿透問題，難以在合理的時間內(nèi)給出可靠的計算結(jié)果。MC分段計算方法［1］雖然在一定程度上提高了MC方法計算效率，但由于分段計算中MC統(tǒng)計誤差的連續(xù)傳遞性，對于有效解決深穿透屏蔽計算問題仍存在困難。

對于先進核裝置及其廠房、反應(yīng)堆堆坑底部、船用堆艙等復(fù)雜屏蔽計算問題，由于其體積龐大，幾何結(jié)構(gòu)復(fù)雜，有大塊屏蔽或復(fù)雜孔道等特點，單一的SN或MC方法無法提供可靠的計算結(jié)果。為了解決這類屏蔽問題，最有效的方法是開發(fā)一種SN方法與MC方法的耦合方法。國際上已經(jīng)進行了這方面的研究，比如美國橡樹嶺實驗室開發(fā)的 DOT3．5-DOMINOMORSE［2］，美 國 加 利 福 尼 亞 大 學(xué) 開 發(fā) 的DORT-PROBGEN-MCNP［3］等。此 類 研 究 工作僅支持二維SN方法與MC方法的耦合，且SN計算只限于R-Z幾何，這無疑難以滿足現(xiàn)代核裝置精確屏蔽設(shè)計的要求。近來日本東芝公司采用DOMINO耦合方法實現(xiàn)了三維離散縱標(biāo)程序 TORT 和 MCNP 的耦合［4］，其 中TORT計算只限于直角坐標(biāo)系，在圓柱坐標(biāo)系問題計算中存在限制。目前國際上還沒用發(fā)布基于三維SN和MC耦合方法的通用程序，而國內(nèi)尚沒有發(fā)現(xiàn)三維SN－MC耦合屏蔽計算方法研究方面的文獻資料。

本文開發(fā)了三維SN－MC耦合計算方法，充分發(fā)揮SN方法解決深穿透問題的優(yōu)勢和MC方法模擬復(fù)雜幾何的長處，提高其屏蔽計算的精度及速度，用于解決大型復(fù)雜核設(shè)施屏蔽計算難題。本文介紹了新開發(fā)的三維耦合屏蔽計算方法、程序系統(tǒng)以及基于該方法的耦合程序系統(tǒng)例題測試。

1 SN－MC耦合方法

SN和MC方法是求解中子輸運問題的兩種不同方法。SN方法是在相空間（r×E×Ω）內(nèi)對空間r、能量E和方向變量Ω采用直接離散方法數(shù)值求解中子輸運方程，得到粒子角注量率ψ（ri，Eg，Ωm）。MC方法通過對大量單個粒子的運動歷史逐個進行跟蹤模擬，然后用統(tǒng)計方法作出隨機變量某個特征的估計量，作為問題的解［5］。一般來說，SN-MC耦合屏蔽計算分析可以遵循下面的步驟：（1）模型劃分：模型分解為適合SN計算的模型區(qū)及適合MC模擬的區(qū)域；（2）指定連接面：定義SN模型和MC模型的連接面；（3）SN計算：采用SN程序?qū)N模型區(qū)進行計算，得到連接面粒子角注量率；（4）轉(zhuǎn)換計算：將SN計算得到的連接面粒子角注量率轉(zhuǎn)換為粒子在空間、能量和方向上的概率分布；（5）MC源粒子抽樣：編寫MC自定義源抽樣程序并嵌入MC程序，根據(jù)轉(zhuǎn)換計算得到的概率分布進行源粒子抽樣，依次確定源粒子的位置、能量和方向；（6）MC計算：采用MC程序，對上述源粒子進行跟蹤模擬計算，得到所需計算結(jié)果。

SN－MC耦合方法關(guān)鍵是SN角注量率到MC源粒子參數(shù)概率轉(zhuǎn)換的實現(xiàn)。三維SN－MC耦合方法采用如下方法實現(xiàn)上述轉(zhuǎn)換：

上式中，ψ（ri，Eg，Ωm）為網(wǎng)格ri、能群 Eg和離散方向Ωm區(qū)間內(nèi)的粒子角注量率；wm為求積權(quán)重系數(shù)；λm表示粒子飛行方向與面法線方向的夾角余弦值；Ai表示網(wǎng)格區(qū)間ri對應(yīng)的連接面面積；I為連接面網(wǎng)格數(shù)目；G為能群數(shù)目；M為離散方向數(shù)目。式（1）p（i）表示粒子落在網(wǎng)格區(qū)間ri內(nèi)的概率；式（2）p（g／i）表示網(wǎng)格區(qū)間ri內(nèi)，粒子落在能群Eg內(nèi)的概率；式（3）p（m／g／i）表示網(wǎng)格區(qū)間ri內(nèi)，能群Eg內(nèi)，粒子落在離散方向Ωm區(qū)間的概率。

2 程序系統(tǒng)

三維SN－MC耦合程序通過編寫接口程序，實現(xiàn)上述耦合方法。耦合程序系統(tǒng)流程圖如圖1所示。程序系統(tǒng)分為4個主要模塊：SN計算模塊、SN-MC接口程序、MC自定義源抽樣程序和MC計算模塊。

圖1 三維SN－MC耦合計算示意圖Fig．1 Sketch of 3D SN－MC coupled calculation

其中三維SN計算得到連接面的粒子角注量率，并以BNDRYS格式存儲輸出。SN-MC接口程序的主要功能是根據(jù)用戶輸入文件，將SN計算輸出的角注量率文件轉(zhuǎn)換為源粒子在空間、能群和方向區(qū)間上的概率分布。用戶輸入文件主要包括以下相關(guān)數(shù)據(jù)：SN模型坐標(biāo)系、連接面網(wǎng)格劃分，能群結(jié)構(gòu)和求積組設(shè)置等。MC-SOURCE源粒子抽樣程序，根據(jù)連接面的粒子概率分布進行MC源粒子參量抽樣，獲得粒子位置、能量和方向信息，為MC計算提供源項。利用MC-SOURCE抽取的源粒子信息進行MC計算，得到復(fù)雜幾何區(qū)的粒子通量分布。該程序系統(tǒng)可以處理三維X-Y-Z及R-θ-Z 幾何。

3 測試與分析

3．1 測試模型描述

采用三維 X-Y-Z 及R-θ-Z 幾何例題對三維SN-MC耦合程序系統(tǒng)進行校核。其中SN計算采用TORT［6］程序及 MATXS10多群截面數(shù)據(jù)庫。MC計算采用MCNP［7］程序及基于ENDF／B-Ⅵ評價核數(shù)據(jù)庫開發(fā)的連續(xù)能量截面數(shù)據(jù)庫［8］。不同于實際應(yīng)用，測試例題選取了單獨采用MCNP程序計算可以獲得精確解的測試模型，以便與耦合計算結(jié)果進行對比分析。兩個測試模型分別如圖2和圖3所示。

圖2 測試模型1Fig．2 Test model 1

模型1是長×寬×高為20 cm×20 cm×24 cm的長方體。20 cm×20 cm×6 cm的中子源位于模型底部，在模型靠近上端水平面沿X方向依次設(shè)置20個探測器。模型2為半徑26 cm，高40 cm的20°圓柱，源區(qū)為半徑6 cm區(qū)域。0°和20°設(shè)為反射邊界條件。其余為真空邊界條件。測試模型的中子源均為各向同性的14 Me V中子，屏蔽材料為不銹鋼和水的混合物（60a%SS316＋40a%Water）。如圖2和圖3所示，連接面將模型分割為兩部分：屏蔽區(qū)和探測器計數(shù)區(qū)。

圖3 測試模型2Fig．3 Test model 2

3．2 結(jié)果比較

耦合計算中，TORT計算建立整個模型，考慮周圍材料對連接面源造成影響的同時，得到連接面角注量率分布和探測器處TORT計算結(jié)果。模型1和模型2分別采用直角坐標(biāo)及圓柱坐標(biāo)建立TORT模型，網(wǎng)格劃分分別為：20×20×24和26×11×21個。TORT計算采用P3階勒讓德展開，S8全對稱高斯求積組，離散方向個數(shù)為96。圖4給出了模型1探測器處MCNP、TORT及TORT-MCNP三種程序計算的中子注量率結(jié)果對比?？梢钥闯?，基于這三種方法的計算結(jié)果吻合得很好。TORTMCNP耦合程序計算結(jié)果總體上介于TORT和 MCNP計算結(jié)果之間。TORT-MCNP與MCNP計算結(jié)果相比，中間探測器結(jié)果偏差在1%以內(nèi)，兩端探測器結(jié)果偏差稍大，在3%左右。這主要是由于SN和MC程序處理邊界條件方式不同，對結(jié)果造成的影響。

圖4 模型1中子注量率比較Fig．4 Comparison of neutron flux of model 1

模型2因為角度方向采用反射邊界條件，同一半徑周向探測器計數(shù)結(jié)果相同，因此選取其中一個探測器進行計算。三種程序計算的探測器中子能譜如圖5所示。TORT-MCNP耦合計算的探測器總中子注量率結(jié)果與MCNP結(jié)果相比相差小于1%?？梢钥闯觯谶@三種方法的計算結(jié)果吻合得很好，顯示三維耦合程序可以提供屏蔽問題的精確解。

圖5 模型2中子能譜比較Fig．5 Comparison of neutron spectrum of model 2

4 結(jié)論

基于SN和MC方法的耦合方法，開發(fā)了TORT-MCNP三維耦合程序系統(tǒng)。該程序集成SN方法解決深穿透問題的優(yōu)點以及MC技術(shù)在復(fù)雜幾何模型模擬方面的長處，克服常規(guī)屏蔽計算方法在大型核裝置屏蔽分析方面的缺點，提高其屏蔽計算的精度及速度。通過X-YZ及R-θ-Z兩種幾何測試例題對耦合程序進行校核計算，并與MCNP及TORT結(jié)果進行比較。結(jié)果吻合良好，初步驗證了三維SN－MC耦合方法及程序系統(tǒng)的正確性。下一步將對三維SN－MC程序系統(tǒng)進行進一步的基準驗證并應(yīng)用于實際工程。

［1］ 鐘兆鵬，施工，胡永明．用MCNP程序計算水平輻照孔道屏蔽 ［J］．清華大學(xué)學(xué)報（自然科學(xué)版），2011，41（12）：16-18．

［2］ Emmett MB，Burgart C E，Hoffman T J．DOMINO，a general purpose code for coupling discrete ordinates and Monte Carlo radiation transport calculations，ORNL-4853［R］．USA：ORNL，1973．

［3］ Eggleston J E，Abdou MA，Youssef MZ．The use of MCNP for neutronics calculations within large buildings of fusion facilities［J］．Fusion Engineering and Design 1998，42：281-288．

［4］ Masahiko Kurosawa，TORT／MCNP coupling method for the calculation of neutron flux around a core of BWR［J］．Radiation Protection Dosimetry，2005，116（1-4）：513-517．

［5］ 謝仲生等．核反應(yīng)堆物理數(shù)值計算［M］．北京：原子能出版社，1997．

［6］ RhoadesW A，Simpson D B．The TORT Three-dimensional Discrete Ordinates Neutron／Photon Trans port Code［R］．ORNL／TM13221，Oak Ridge National Laboratory，1997．

［7］ BRIESMEISTER J F．MCNP：A general Monte Carlo N-particle transport code，version 4C，LA-13709-M［R］．USA：LANL，2000．

［8］ WIENKE H，HERMAN M．FENDL／MG-2．0 and FENDL／MC-2．0，the processed cross-section libraries for neutronphoton transport calculations，version 1［R］．Vienna，Austria：Nuclear Data Section，International Atomic Energy Agency，1998．