基于離散縱標法與蒙特卡羅方法的三維耦合程序開發

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

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

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

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

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

1 SN－MC耦合方法

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

SN－MC耦合方法關鍵是SN角注量率到MC源粒子參數概率轉換的實現。三維SN－MC耦合方法采用如下方法實現上述轉換：

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

2 程序系統

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

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

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

3 測試與分析

3．1 測試模型描述

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

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

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

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

3．2 結果比較

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

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

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

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

4 結論

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

［1］ 鐘兆鵬，施工，胡永明．用MCNP程序計算水平輻照孔道屏蔽 ［J］．清華大學學報（自然科學版），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］ 謝仲生等．核反應堆物理數值計算［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．