基于蒙特卡羅-離散縱標雙向耦合方法的快中子注量基準分析

袁龍軍，陳義學，2，韓靜茹

（1.華北電力大學 核科學與工程學院，北京 102206；2.國家核電技術公司 北京軟件技術中心，北京 100029；3.環境保護部 核與輻射安全中心，北京 100082）

蒙特卡羅法（MC）和離散縱標法（SN）是最為常用的屏蔽計算方法。蒙特卡羅法可精確模擬復雜幾何模型，但計算耗時，尤其對于屏蔽計算中常見的深穿透問題。與蒙特卡羅法相比，離散縱標法計算速度快，但難以精確描述復雜的幾何模型。為綜合兩種方法的優點以解決同時具有復雜幾何模型和深穿透問題的屏蔽計算問題，國內外開展了基于蒙特卡羅法和離散縱標法的耦合方法的研究，研制了大量耦合計算程序。目前國內外公開發表的一維耦合程序有HETC96-ANISN［1］、HERMES-ANISN［2］、MCNPANSIN［3］等；二 維 耦 合 程 序 有MCNP-DORT［3］、DORT-MCNP［4］、DORT-PROBEGEN-MCNP［5］、DORT-DOMINO-MORSE［6］等；三維耦合程序有MCNP-TRIDENT［7］、TORT-MCNP［8］等。然而一維和二維耦合程序不能精確計算三維模型，而三維耦合程序僅為單向耦合，在工程計算中不能根據實際情況靈活使用蒙特卡羅法和離散縱標法，嚴重限制了耦合方法在實際工程中的應用。

為解決上述問題，開發基于蒙特卡羅法和離散縱標法的三維蒙特卡羅-離散縱標雙向耦合計算方法，將耦合方法用于壓水堆壓力容器快中子注量計算，并與基準報告進行比較分析。

1 MC-SN 雙向耦合方法

作為輻射屏蔽計算最為常用的兩種計算方法，蒙特卡羅法和離散縱標法分別為非確定論方法和確定論方法。實現蒙特卡羅法和離散縱標法耦合計算的關鍵在于實現蒙特卡羅法模擬的粒子徑跡信息和離散縱標法計算的角通量密度間的相互轉換，即實現蒙特卡羅模擬粒子的位置、能量和飛行方向與離散縱標法計算的角通量密度的空間網格、能群和離散方向的一一對應轉換。

其中，MC粒子徑跡信息到SN角通量密度的轉換關系［9］為：

式中：ψi，j，k，m，g為空 間 網 格（i，j，k）處 第g 群 內的m 方向范圍內的中子角通量密度；weightn為MC 粒子n 的權重；wm為求積權重系數；N 為MC模擬的源粒子數；ΔA 為給定面元的面積；λn為MC粒子徑跡和面法線方向夾角的余弦值；En、rn和Ωn分別為粒子的能量、空間位置和飛行方向。

SN角通量密度到MC 粒子徑跡信息的轉換關系［10］為：

式中：wm為求積權重系數；ψ（ri，Eg，Ωm）為網格ri、能群Eg和離散方向Ωm區間內的粒子角通量密度；λm為粒子飛行方向與面法線方向夾角的余弦值；Ai為網格區間ri對應的連接面的面積；求和上限I、G 和M 分別為連接面網格數、能群數和離散方向數；p（ri）為MC抽樣粒子位于網格區間ri內的概率；p（Eg／ri）為在網格區間ri內抽樣、粒子位于能群區間Eg內的概率；p（Ωm／Eg／ri）為在網格區間ri、能群區間Eg內抽樣，粒子位于離散方向Ωm內的概率。

根據以上轉換方法，通過自主研發的接口程序，實現蒙特卡羅法和離散縱標法的雙向耦合，其流程圖示于圖1，其中MC 程序采用國際通用的粒子輸運程序MCNP4C［11］，并采用ENDF60連續能量截面數據庫；SN程序采用美國橡樹嶺國家實驗室開發的三維離散縱標程序TORT［12］，采用多群截面庫MATXS10。

2 基準驗證

為驗證三維MC-SN雙向耦合屏蔽計算方法的正確性和可靠性，本文采用了BNL（Brookhaven National Laboratory）開 發 的 NUREG／CR-6115［13］壓水堆壓力容器快中子注量計算基準例題。

基準題的壓水堆模型如圖2所示。反應堆采用基準題中標準的IN-OUT 堆芯裝載模式；堆芯采用204 個燃耗深度不同的組件?；贛C-SN雙向耦合方法，并根據基準題模型的特點，選取熱屏蔽內表面作為第1耦合面、壓力容器堆焊層作為第2耦合面，將模型劃分為MC模擬區和SN模擬區，如圖2a所示。其計算流程如下。

圖1 三維蒙特卡羅-離散縱標雙向耦合方法流程圖Fig.1 Flow chart of program system for 3D MC-SN bidirectional coupling method

圖2 1／8壓水堆模型水平（a）和垂直（b）剖面圖Fig.2 Horizontal（a）and vertical（b）sections of 1／8PWR model

1）MCNP計算。利用MCNP程序實現堆芯到熱屏蔽內表面區域的精確計算，得到第1耦合面的中子徑跡信息。

2）MCNP-TORT 計算。通過接口程序，將MCNP程序的粒子徑跡信息轉換為TORT程序的中子角通量密度，得到用于TORT 計算的邊界源文件。

3）TORT 計算。通過邊界源文件，利用TORT 程序實現熱屏蔽內表面至壓力容器焊層的計算。

4）TORT-MCNP 計算。通過接口程序，將TORT 程序的中子角通量密度轉換為MCNP程序的粒子概率分布，得到用于MCNP計算的粒子抽樣文件。

5）MCNP計算。通過粒子抽樣文件，利用MCNP程序實現壓力容器內部計算區域的精確計算。

通過以上流程，實現MCNP-TORT 雙向耦合計算。

根據基準題模型結構的對稱性，本文選取了全堆的1／8作為計算模型，在0°和45°邊界采用反射邊界，在壓力容器外表面和模型頂部及底部采用真空邊界。其中壓力容器內壁峰值、壓力容器1／4壁厚峰值和壓力容器內壁焊縫處為計數位置，如圖2b所示。

在基準報告中，選取壓力容器內壁峰值、1／4壁厚峰值和內壁焊縫處的快中子注量率（E＞1.0 MeV）作為計算區域，并給出MCNP4A及DORT 的計算結果。本文使用單獨的MCNP4C 程序和MCNP-TORT 雙向耦合程序進行相應計算，得到了對應的計算結果，并與基準報告進行對比，如圖3所示。

從圖3可看出，耦合程序計算結果與基準結果符合良好，平均相對誤差不超過3%；與MCNP4A結果相比，最大相對誤差不超過10%。

通過對比，較SN程序，MCNP-TORT 雙向耦合程序計算結果與基準結果更為吻合。

在相同的計算時間下，MCNP-TORT 雙向耦合程序較單一MCNP程序統計誤差更小，計算結果可靠性更強。單一MCNP 程序計算耗時53min，模擬粒子數為1.5 億，統計誤差為3%左右。MCNP-TORT 雙向耦合程序的計算時間為48min，其中兩次MCNP計算分別耗時12min和23min，模擬粒子數分別為5 000萬和9 000萬，統計誤差均為1%左右；TORT 計算耗時13min。

圖3 壓力容器內壁峰值、1／4壁厚峰值和內壁焊縫處的快中子注量率（E＞1.0 MeV）周向分布Fig.3 Circular distributions of E＞1.0 MeV fast neutron fluence rate at pressure vessel inner wall axial peak location，pressure vessel 1／4axial peak location and pressure vessel lower weld

3 結論

為解決同時具有復雜幾何和深穿透問題的屏蔽計算問題，開發了基于三維蒙特卡羅-離散縱標雙向耦合方法。本文將此耦合計算方法應用于NUREG／CR-6115壓水堆壓力容器快中子注量計算，并與基準報告中其他計算方法進行對比，結果符合良好；在相同的計算時間下，耦合計算程序統計誤差優于單一蒙特卡羅程序，可靠性更強，驗證了三維蒙特卡羅-離散縱標雙向耦合方法的正確性與可靠性。

［1］ GABRIEL T A.CALOR：A Monte Carlo program package for the design and analysis of calorimeter systems，ORNL／TM-5619［R］.USA：ORNL，1977.

［2］ CLOTH P，FILGES D，NEEF R D，et al.HERMES：A Monte Carlo program system for beam-materials interaction studies［R］.Germany：KFA，1988.

［3］ GALLMEIER F X，PEVEY R E.Creation of a set of interface utilities to allow coupled Monte Carlo／discrete ordinate shielding analysis［C］∥Third International Topical Meeting on Nuclear Applications of Accelerator Technology.USA：American Nuclear Society，1999：404-409.

［4］ 鄭征，吳宏春，曹良志.蒙卡與SN耦合方法在壓水堆堆腔漏束計算中的應用［C］∥第十一屆全國蒙特卡羅方法及應用學術交流會.綿陽：［出版者不詳］，2012.

［5］ 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：281-288.

［6］ 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：Oak Ridge National Laboratory，1973.

［7］ URBAN W T，SEED T J，DUDZIAK D J.Nucleonic analysis of a preliminary design for the ETF neutral-beam-injector duct shielding，LAUR-80-2926［R］.New Mexico：Los Alamos Scientific Laboratory，1980.

［8］ 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.

［9］ CHEN Y，FISCHER U.Program system for three-dimensional coupled Monte Carlo-deterministic shielding analysis with application to the accelerator-based IFMIF neutron source［J］.Nuclear Instruments and Methods in Physics Research A，2005，551：387-395.

［10］HAN Jingru，CHNE Yixue，YUAN Longjun，et al.Validation of three-dimensional coupled discrete ordinates-Monte Carlo program system through shielding benchmark analysis［R］.［S.l.］：ICETCE，2012.

［11］BRIESMEISTER J F.MCNP：A general Monte Carlo N-particle transport code，version 4C，LA-13709-M［R］.USA：Los Alamos National Laboratory，2000.

［12］RHOADES W A，SIMPSON D B.The TORT three-dimensional discrete ordinates neutron／photon transport code，ORNL／TM13221［R］.USA：Oak Ridge National Laboratory，1997.

［13］CAREW J F，HU K，AROSON A，et al.PWR and BWR pressure vessel fluence calculation benchmark problem and solutions［R］.Washington D.C.：Brookhaven National Laboratory，2001.