三維離散縱標-蒙特卡羅耦合系統TDOMINO開發與驗證

韓靜茹，陳義學，袁龍軍，張春明

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

在求解輻射屏蔽問題時，常用的蒙特卡羅方法和離散縱標法各有長處和局限性。三維離散縱標(SN)-蒙特卡羅(MC)耦合是利用粒子角通量密度分布和MC源變量抽樣概率之間的轉換算法，將三維離散縱標程序與蒙特卡羅程序耦合，充分發揮各自的優勢，同時又克服各自的局限性。對大型復雜屏蔽問題而言，三維耦合計算可更加有效地解決深穿透問題和詳細描述復雜幾何，揭示更為真實的安全裕量，有利于核電站安全性的改進和經濟性的提高。自20世紀70年代至今，國內外科研人員對離散縱標-蒙特卡羅耦合計算進行了大量的研究，如美國橡樹嶺國家實驗室開發的DOMINO[1]、美國加利福尼亞大學開發的PROBGEN[2]、國內西安交通大學開發的DORT2MCNP[3]等，但上述研究僅支持r-z幾何下的二維離散縱標計算和蒙特卡羅計算的耦合。日本東芝公司實現了三維離散縱標程序TORT和MCNP的耦合[4]，但僅限于x-y-z幾何，未解決r-θ-z幾何問題。這嚴重限制了耦合方法在實際工程中的應用。因此，本文進行三維離散縱標-蒙特卡羅耦合系統TDOMINO的開發，并針對國際通用基準題進行驗證計算。

1 三維SN-MC耦合系統框架

TDOMINO耦合系統流程圖示于圖1。以現有程序為基礎，通過開發建立相應的耦合接口和源粒子抽樣程序來實現程序數據交換和計算。這樣的耦合框架既避免了對現有程序結構的修改，也便于今后對程序進行維護和升級。同時也使耦合形式具有較好的靈活性，可根據問題分析的需要選擇用于耦合計算的程序。構成TDOMINO耦合系統的基礎程序包括：三維離散縱標程序TORT[5]和蒙特卡羅程序MCNP[6]。在耦合系統中，為充分發揮兩種程序的優勢，TORT和MCNP分別用以模擬大塊屏蔽區和復雜幾何區計算。其中三維SN計算得到連接面的粒子角通量密度，并以BNDRYS格式存儲輸出，作為SN-MC接口程序的輸入文件之一。SN-MC接口程序將連接面的角通量密度轉換為粒子概率分布，再經MC-SOURCE源粒子抽樣程序轉換為MC源粒子信息參數，為MC計算提供源項，實現三維SN和MC耦合輸運計算，得到復雜幾何區的粒子通量分布。該程序系統可處理三維x-y-z及r-θ-z幾何。

圖1 TDOMINO耦合系統流程圖

2 三維SN-MC耦合方法

對于上述兩個基礎程序，由于其物理建模和數值求解方法不同，耦合需解決連接區內SN計算的粒子角通量密度分布和MC計算所需源項之間映射的問題。在TDOMINO耦合系統中，通過開發三維SN-MC接口程序和MC-SOURCE源粒子抽樣程序分別實現連接區內SN計算得到的粒子角通量密度和粒子概率分布之間映射及粒子概率分布和MC源粒子信息間的轉換。具體算法為：

(1)

(2)

(3)

(4)

式中：ψ(ri，Eg，Ωm)為與空間、能量和方向相關的粒子角通量密度；i為空間網格(總網格數用I表示)；g為能群 (總能群數用G表示)；m為離散方向(總離散方向數用M表示)；Ai為第i個空間網格的面積；λm為SN求積組的方向余弦μm、ξm或ηm；wm為λm對應的SN求積組的權重；p(i)、p(g/i)和p(m/g/i)分別為粒子落在空間網格區間ri內的概率、網格區間ri內粒子能量落在能群Eg內的概率和網格區間ri、能群Eg內粒子落在離散方向m內的概率。

MC-SOURCE源粒子抽樣程序基于三維SN-MC接口程序計算得到的粒子概率分布，根據式(4)進行隨機抽樣，依次獲得粒子所在的SN計算劃分的網格Δr、能群ΔE和離散方向ΔΩ區間，并在這些子區內均勻抽樣，最終依次確定源粒子的位置(x,y,z)、能量(E)及飛行方向與(x,y,z)3個坐標軸的夾角余弦值(u，v，w)等信息，為下一步MC計算提供源項。需注意，圓柱坐標系下的源粒子位置和飛行方向表達方式與直角坐標系下的不同，需進行相應的坐標轉換。

3 耦合程序的驗證

為驗證三維SN-MC耦合程序，采用美國橡樹嶺國家實驗室公布的HBR-2屏蔽基準題[7]，利用TDOMINO耦合系統進行驗證計算，并針對圓柱坐標系和直角坐標系下的耦合模型特點，采用不同的耦合形式進行計算。

三維SN-MC耦合計算可分為以下步驟。1) 模型劃分：分為適合SN計算的屏蔽區及適合MC模擬的復雜幾何區；2) 指定公共幾何:連接大塊屏蔽區(SN網格區)和復雜幾何區(MC模型)的公共面；3) SN計算：得到公共面內每個SN網格的角粒子通量密度分布；4) 轉換計算：將SN網格的角通量密度分別轉換為粒子在空間、能量和方向上的分布；5) MC源抽樣概率：將粒子分布轉換成MC源參數相應的概率分布，并嵌入MC源抽樣程序；6) MC計算：根據上步得到的源參數分布概率進行抽樣模擬，進行MC計算。

3.1 HBR-2基準題

HBR-2基準題是在美國橡樹嶺國家實驗室發布的用于驗證壓水堆壓力容器屏蔽計算例題。HBR-2是熱功率為2 300 MW的壓水堆核電機組，反應堆堆芯包含157個燃料組件，壓力容器內半徑為197.485 cm。堆芯外依次圍繞有堆芯圍板、吊籃、熱屏蔽、輻照監督管、壓力容器以及生物屏蔽等部件，更詳細的描述見HBR-2基準報告[7]。輻照監督管中心與堆芯中心的距離為191.15 cm，設置在熱屏蔽外壁，位于20°方位角處。基準實驗在反應堆第9循環內輻照監督管處放置測量儀對試樣比活度進行測量。基準報告中提供了實驗測量值和1/8模型的DORT程序計算結果。

3.2 計算模型

根據基準題提供的具體參數及反應堆對稱性，結合TDOMINO應用特點，分別建立r-θ-z幾何模型和x-y-z幾何模型進行耦合計算。堆芯部分采用三維SN程序TORT建模計算，熱屏蔽外的輻照監督管采用MCNP程序建立精確模型。

r-θ-z坐標下三維SN-MC耦合計算模型示于圖2。以基準報告中計算模型為例，建立1/8堆型模型，0°和45°的角邊界設為反射邊界，壓力容器外表面設為真空邊界，同時將堆芯上反射層和下反射層的外表面設為軸向真空邊界條件。將吊籃外表面設為SN-MC耦合模型的連接面，考慮到中子的散射作用，SN模型和MC模型有一段重疊區。堆芯到熱屏蔽設為SN模型，采用TORT圓柱坐標系建模計算，吊籃外表面到壓力容器設為MC模型，采用MCNP程序詳細描述輻照監督管模型及材料，并在輻照監督管處設置了中子注量率計數卡。TDOMINO計算中，TORT和MCNP程序分別采用基于ENDF/B-Ⅵ的TEXT-10多群截面庫和點截面庫。耦合計算中TORT在r、θ、z三個方向上的網格個數分別取為118、52、88。

圖2 三維SN-MC耦合計算r-θ-z模型(1/8堆芯)

為驗證TDOMINO在x-y-z幾何中的應用及處理多面耦合的能力，結合TORT建模特點，建立的1/4 HBR-2全堆芯組件均勻模型為三維SN-MC耦合計算模型，如圖3所示。將圍板設為耦合面，堆芯到圍板的SN模型區采用TORT程序建立x-y-z幾何模型(圖3b)，圍板到壓力容器采用MCNP程序建模計算(圖3c)，實現直角坐標系下的三維SN-MC多面耦合計算。

a——耦合計算模型；b——TORT計算模型；c——MCNP計算模型

3.3 計算結果

分別采用三維SN-MC耦合程序TDOMINO和單獨的MCNP程序對HBR-2基準模型輻照監督管處6種典型核素的比活度進行計算，并與實驗測量、離散縱標法程序計算結果進行對比分析。

表1列出HBR-2基準模型輻照監督管處6種典型核素比活度的實驗測量值。表2列出不同程序計算的6種核素比活度與實驗測量值間的比值(C/M)。其中實驗測量值和DORT采用BUGLE-96[8]庫的計算結果取自HBR-2基準報告，TORT結合TEXT10多群截面數據庫的計算結果取自文獻[9]，SN-MC和SN-MC-2分別表示TDOMINO在圓柱坐標系和直角坐標系下的計算結果。MCNP表示堆芯組件均勻模型的單一MC計算結果，MCNP-PIN表示堆芯PIN功率的計算結果取自文獻[10]。DORT、SN-MC、SN-MC-2、MCNP、MCNP-PIN和TORT程序計算得到的輻照監督管處的平均C/M值分別為0.89±0.04、1.03±0.04、1.04±0.03、1.10±0.04、0.91±0.03和1.04±0.04。計算結果表明，TDOMINO耦合系統計算結果與實驗測量值和其他程序計算值吻合較好，初步驗證了三維SN-MC耦合系統TDOMINO在不同空間坐標系中應用的正確性。在上述計算中，SN-MC耦合程序的計算時間為3 908 min，其中，SN-MC中MCNP計算耗時3 869 min，最大統計誤差為3%左右，SN-MC中TORT計算耗時39 min。單一MCNP程序計算耗時7 180 min，最大統計誤差為5%左右。表明三維SN-MC耦合程序與單一MCNP程序相比，在較少的時間下得到的統計誤差更小，計算結果可靠性更強。

表1 實驗測量所得的比活度

表2 計算與測量的比活度之比

4 結論

為校核驗證三維離散縱標-蒙特卡羅耦合程序系統在復雜模型輻射屏蔽計算問題中的應用可行性，針對美國HBR-2壓力容器基準實驗，對輻照監督管處的中子能譜和幾種典型核素的比活度進行了計算分析。采用三維TDOMINO耦合系統分別建立了圓柱坐標系下的HBR-2 1/8模型和直角坐標系下的HBR-2 1/4耦合計算模型，并進行了驗證計算。在輻照監督管處計算的平均C/M值分別為1.03±0.04和1.04±0.03，可看出耦合計算結果與實驗測量及其他程序計算值吻合良好，證明了TDOMINO耦合系統開發的正確性。下一步將針對具體的工程實際問題(如反應堆堆腔漏束輻射計算等)，開展三維離散縱標-蒙特卡羅耦合系統的應用研究。

參考文獻：

[1] EMMETT M B, 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.

[2] EGGLESTON J E, ABDOU M A, YOUSSEF M Z. The use of MCNP for neutronics calculations within large buildings of fusion facilities[J]. Fusion Engineering and Design, 1998, 42(1-4): 281-288.

[3] 鄭征，吳宏春，曹良志，等. 蒙特卡羅與離散縱標耦合方法在壓水堆堆腔漏束計算中的應用[J]. 強激光與粒子束,2012,24(12):2 946-2 950.

ZHENG Zheng, WU Hongchun, CAO Liangzhi,et al. Application of Monte Carlo and discrete ordinate coupling method to pressurized water reactor cavity radiation streaming calculation[J]. High Power Laser and Particle Beams, 2012, 24(12): 2 946-2 950(in Chinese).

[4] KUROSAWA M. 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] BRIESMEISTER J F. MCNP: A general Monte CarloN-particle transport code, Version 4C, LA-13709-M[R]. USA: LANL, 2000.

[6] RHOADES W A, SIMPSON D B. The TORT three-dimensional discrete ordinates neutron/photon transport code, ORNL/TM13221[R]. USA: ORNL, 1997.

[7] REMEC F B K, KAM H B. Robinson-2 pressure vessel benchmark[R]. USA: ORNL, 1997.

[8] WHITE J E. BUGLE-96: Coupled 47 neutron, 20 gamma-ray group cross section library derived from ENMF/B-Ⅵ for LWR shielding and pressure vessel dosimetry applications[R]. USA: RSIC Data Library Collection, 1996.

[9] 楊壽海,陳義學,王偉金,等. 三維離散縱標方法在RPV快中子注量率計算中的初步應用[J]. 核科學與工程，2011，31(4)：294-298.

YANG Shouhai, CHEN Yixue, WANG Weijin, et al. The analysis of RPV fast neutron flux calculation for PWR with three-dimensional SNmethod[J]. Chin J Nucl Sci Eng, 2011, 31(4): 294-298(in Chinese).

[10] 石生春. 基于蒙特卡羅方法的壓水堆壓力容器快中子注量率的計算分析[D]. 北京：華北電力大學，2010.