基于離散縱標(biāo)法與蒙特卡羅方法的三維耦合程序開(kāi)發(fā)

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

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

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

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

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

1 SN－MC耦合方法

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

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

上式中，ψ（ri，Eg，Ωm）為網(wǎng)格ri、能群 Eg和離散方向Ωm區(qū)間內(nèi)的粒子角注量率；wm為求積權(quán)重系數(shù)；λm表示粒子飛行方向與面法線(xiàn)方向的夾角余弦值；Ai表示網(wǎng)格區(qū)間ri對(duì)應(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耦合程序通過(guò)編寫(xiě)接口程序，實(shí)現(xiàn)上述耦合方法。耦合程序系統(tǒng)流程圖如圖1所示。程序系統(tǒng)分為4個(gè)主要模塊：SN計(jì)算模塊、SN-MC接口程序、MC自定義源抽樣程序和MC計(jì)算模塊。

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

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

3 測(cè)試與分析

3．1 測(cè)試模型描述

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

圖2 測(cè)試模型1Fig．2 Test model 1

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

圖3 測(cè)試模型2Fig．3 Test model 2

3．2 結(jié)果比較

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

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

模型2因?yàn)榻嵌确较虿捎梅瓷溥吔鐥l件，同一半徑周向探測(cè)器計(jì)數(shù)結(jié)果相同，因此選取其中一個(gè)探測(cè)器進(jìn)行計(jì)算。三種程序計(jì)算的探測(cè)器中子能譜如圖5所示。TORT-MCNP耦合計(jì)算的探測(cè)器總中子注量率結(jié)果與MCNP結(jié)果相比相差小于1%?？梢钥闯觯谶@三種方法的計(jì)算結(jié)果吻合得很好，顯示三維耦合程序可以提供屏蔽問(wèn)題的精確解。

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

4 結(jié)論

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

［1］ 鐘兆鵬，施工，胡永明．用MCNP程序計(jì)算水平輻照孔道屏蔽 ［J］．清華大學(xué)學(xué)報(bào)（自然科學(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ù)值計(jì)算［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．