袁龍軍，陳義學(xué)，2，韓靜茹

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

蒙特卡羅法（MC）和離散縱標(biāo)法（SN）是最為常用的屏蔽計算方法。蒙特卡羅法可精確模擬復(fù)雜幾何模型，但計算耗時，尤其對于屏蔽計算中常見的深穿透問題。與蒙特卡羅法相比，離散縱標(biāo)法計算速度快，但難以精確描述復(fù)雜的幾何模型。為綜合兩種方法的優(yōu)點(diǎn)以解決同時具有復(fù)雜幾何模型和深穿透問題的屏蔽計算問題，國內(nèi)外開展了基于蒙特卡羅法和離散縱標(biāo)法的耦合方法的研究，研制了大量耦合計算程序。目前國內(nèi)外公開發(fā)表的一維耦合程序有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］等。然而一維和二維耦合程序不能精確計算三維模型，而三維耦合程序僅為單向耦合，在工程計算中不能根據(jù)實際情況靈活使用蒙特卡羅法和離散縱標(biāo)法，嚴(yán)重限制了耦合方法在實際工程中的應(yīng)用。

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

1 MC-SN 雙向耦合方法

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

其中，MC粒子徑跡信息到SN角通量密度的轉(zhuǎn)換關(guān)系［9］為：

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

SN角通量密度到MC 粒子徑跡信息的轉(zhuǎn)換關(guān)系［10］為：

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

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

2 基準(zhǔn)驗證

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

基準(zhǔn)題的壓水堆模型如圖2所示。反應(yīng)堆采用基準(zhǔn)題中標(biāo)準(zhǔn)的IN-OUT 堆芯裝載模式；堆芯采用204 個燃耗深度不同的組件。基于MC-SN雙向耦合方法，并根據(jù)基準(zhǔn)題模型的特點(diǎn)，選取熱屏蔽內(nèi)表面作為第1耦合面、壓力容器堆焊層作為第2耦合面，將模型劃分為MC模擬區(qū)和SN模擬區(qū)，如圖2a所示。其計算流程如下。

圖1 三維蒙特卡羅-離散縱標(biāo)雙向耦合方法流程圖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程序?qū)崿F(xiàn)堆芯到熱屏蔽內(nèi)表面區(qū)域的精確計算，得到第1耦合面的中子徑跡信息。

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

3）TORT 計算。通過邊界源文件，利用TORT 程序?qū)崿F(xiàn)熱屏蔽內(nèi)表面至壓力容器焊層的計算。

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

5）MCNP計算。通過粒子抽樣文件，利用MCNP程序?qū)崿F(xiàn)壓力容器內(nèi)部計算區(qū)域的精確計算。

通過以上流程，實現(xiàn)MCNP-TORT 雙向耦合計算。

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

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

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

通過對比，較SN程序，MCNP-TORT 雙向耦合程序計算結(jié)果與基準(zhǔn)結(jié)果更為吻合。

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

圖3 壓力容器內(nèi)壁峰值、1／4壁厚峰值和內(nèi)壁焊縫處的快中子注量率（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 結(jié)論

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

［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，F(xiàn)ILGES 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耦合方法在壓水堆堆腔漏束計算中的應(yīng)用［C］∥第十一屆全國蒙特卡羅方法及應(yīng)用學(xué)術(shù)交流會.綿陽：［出版者不詳］，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，F(xiàn)ISCHER 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.