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

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

屏蔽系統(tǒng)設計是核裝置工程設計的核心內容之一，其設計優(yōu)劣直接影響工程造價及工作人員與周圍環(huán)境的輻射安全。選取合適的屏蔽計算方法是保證屏蔽系統(tǒng)設計質量的關鍵。華北電力大學在輻射屏蔽方法及應用方面開展了大量的研究工作[1－5］。蒙特卡羅方法（MC）和離散縱標方法（SN）是最常用的屏蔽計算方法。MC方法的優(yōu)點是可以精確模擬復雜幾何模型，但計算耗時，尤其對于屏蔽計算中常見的深穿透問題，有的甚至根本無法得出計算結果。盡管有MC分段計算[6］方法的研究用于解決深穿透的問題，但由于分段計算中MC統(tǒng)計誤差的連續(xù)傳遞性，但這仍舊是其在研究和實際應用中最大的困難。與MC方法相比，SN方法特別適合解決深穿透問題但卻難以精確描述復雜幾何模型。對于大型復雜核裝置屏蔽計算問題，MC方法或SN無法提供可靠的計算結果。

為了解決這類具有復雜物理過程與幾何結構，同時具有深穿透特點的屏蔽問題，結合MC和SN方法優(yōu)勢的耦合計算方法成為首選。如基于MC和SN方法開發(fā)的耦合程序MCNPANSIN[7］、 HETC96－ANISN[8］、 HERMESANISN[9］、MCNP－DORT[7］、MCNP－TRIDENT[10］等。然而這些程序僅支持MC方法與一維或二維離散縱標法的耦合，不支持三維屏蔽計算，難以滿足現(xiàn)代核裝置精確屏蔽設計的要求。陳義學等開發(fā)了三維 MC－SN耦合方法，并用于直角坐標系下國際聚變材料輻照裝置（IFMIF）屏蔽設計計算[11－12］。但該耦合方法在圓柱坐標系中應用時存在限制。本工作對三維MC－SN耦合方法進行改進，然后將其用于圓柱坐標系下PWR壓力容器快中子注量計算，并與基準報告結果進行比較分析，驗證三維 MC－SN耦合方法計算結果的可靠性，從而保證為屏蔽系統(tǒng)優(yōu)化設計的質量提供有力的技術支持。

1 MC－SN耦合方法及程序

MC和SN方法是求解中子輸運問題的兩種不同方法。其中MC方法屬于非確定論方法，通過對大量中子行為觀察分析，用統(tǒng)計平均的辦法，推測出估計量之解。而離散縱標SN方法屬于確定論方法，是用數(shù)值方法求解玻爾茲曼輸運方程的重要方法之一。為了精確處理大型復雜核裝置屏蔽問題，結合MC和SN方法的優(yōu)勢，開發(fā)了三維MC－SN耦合方法。一般來說，耦合屏蔽計算分析可以遵循下面的步驟：（1）模型分解為適合 Monte Carlo模擬的源區(qū)或復雜幾何模型及屏蔽區(qū)SN模型；（2）定義Monte Carlo模型與SN模型的連接面；（3）對源模型進行Monte Carlo模擬，得到通過連接面的粒子徑跡信息；（4）將記錄的Monte Carlo粒子徑跡轉化為SN角通量分布；（5）根據(jù)得到的角通量分布生成SN程序所需的面源文件；（6）利用面源進行SN計算，得到厚屏蔽區(qū)的粒子通量分布。

為了耦合兩種不同的方法，關鍵是MC模擬的公共面粒子徑跡到SN角通量密度的轉換，即將MC粒子的位置、能量和飛行方向分別與SN空間網(wǎng)格、能群和離散方向一一對應轉換。三維 MC－SN耦合方法采用映射方法[11－12］實現(xiàn)上述轉換。已有的三維 MC－SN 耦合方法只考慮了直角坐標系下MC粒子飛行方向與SN離散方向的轉換，在圓柱坐標系下卻不適用。本工作針對上述問題進行了如下改進。MC描述幾何時以直角坐標系為基礎，MC粒子徑跡包括粒子飛行方向分別與x，y，z坐標軸的夾角余弦值u，v，w，通過三個夾角余弦值可以確定粒子運動方向在直角坐標系中對應的極角θ和方位角φ。而在圓柱坐標系中，SN離散方向區(qū)域可由極角θ和ω來表示。其中ω是粒子方向和Z軸形成的平面與R和θ形成的平面間的夾角，可認為是圓柱坐標系中對應的方位角。當點的位置及其空間坐標R改變時，R和θ的坐標軸方向亦隨之發(fā)生改變。將圓柱坐標系中的方位角ω轉換為直角坐標系中的方位角φ，實現(xiàn)SN圓柱坐標系下離散方向與MC徑跡中的粒子方向對應轉換。同時本工作將原三維MC－SN耦合方法中圓柱坐標系幾何范圍由（－90°，90°）擴展到（0°，360°）。

三維MC－SN耦合計算程序系統(tǒng)流程圖如圖1所示。其中Monte Carlo模擬使用國際通用的粒子輸運程序MCNP[13］，三維離散縱標法計算使用美國橡樹嶺國家實驗室開發(fā)的三維離散縱標法程序TORT[14］，而接口程序則基于上述介紹的三維映射方法，將MCNP計算得到的粒子徑跡轉換為TORT計算所需的邊界源文件。

圖1 MC－SN三維耦合計算程序系統(tǒng)流程圖Fig.1 Flow chart of the program system for three－dimensional coupled MC－SN

2 MC－SN基準測試

三維MC－SN耦合屏蔽計算方法采用BNL（Brookhaven National Laboratory）開 發(fā) 的NUREG／CR－6115[15］壓水堆基準例題進行驗證計算，并與基準報告提供的 MCNP、DORT結果進行比較分析。

2.1 基準模型描述

基準計算對象選取 NUREG／CR－6115（BNLNUREG－52395）中標準的IN－OUT 堆芯裝載模式的反應堆。其堆芯由204個燃耗深度不同的組件組成，堆芯外依次由圍板、吊籃、熱屏蔽、壓力容器等結構組成。堆芯高度為335.28cm，堆芯上底面反射層厚度為32.865cm，堆芯下底面反射層厚度為13.97cm。壓力容器內半徑為219.075cm，厚21.59cm，包括0.635cm內表面不銹鋼覆蓋層。

選取熱屏蔽內表面作為耦合計算的公共交界面，將模型劃分為MC模擬區(qū)和SN模擬區(qū)。堆芯到熱屏蔽使用MCNP程序計算，為了減少外圍組件不可忽略的徑向功率梯度對結果的影響，建立堆芯外圍3層組件15×15的pin－bypin精細模型。同時考慮粒子散射對交界面源造成的影響，將MC模型擴建到壓力容器內表面。熱屏蔽到壓力容器部分屬于模型相對簡單的SN模擬區(qū)，采用TORT進行建模計算。接口程序的主要功能是將MCNP計算獲得穿過熱屏蔽內表面的中子徑跡信息（包括權重、空間位置、能量、飛行方向等）轉換為TORT邊界源文件，用于SN模擬區(qū)屏蔽計算，實現(xiàn)三維MCSN耦合計算。

選取與基準報告中相同的計算模型，即1／8反應堆模型。0°和45°處采用反射邊界，壓力容器外表面設置真空邊界，同樣在頂部和底部反射層的外表面采用軸向真空邊界條件。圖2a和2b分別給出了三維耦合計算幾何模型的水平和垂直剖面圖。壓力容器內壁峰值位置為R＝219.393cm，Z＝125.488cm；壓力容器1／4峰值位置為R＝224.473cm，Z＝125.488cm；壓力容器內壁焊縫位置為R＝219.393cm，Z＝67.104 8cm，具體如圖2b所示。

圖2 a 1／8壓水堆模型水平剖面圖Fig.2a Horizontal section of 1／8PWR model

2.2 計算結果及分析

MCNP計算準確地描述計算問題的幾何結構及材料，采用ENDF60截面數(shù)據(jù)庫，它是基于ENDF／B－Ⅵ評價核數(shù)據(jù)庫開發(fā)的連續(xù)能量截面數(shù)據(jù)庫。堆芯內區(qū)組件采用組件平均功率，堆芯外圍組件15×15的精細模型，并采用pin－by－pin的精細功率分布。堆芯瞬發(fā)中子采用混合裂變譜，考慮燃耗對多種裂變同位素的影響。TORT計算采用R－θ－Z幾何模型，P3階勒讓德展開和S8全對稱高斯求積組。R、θ、Z三個方向分別劃分了45、20和38個網(wǎng)格。同時計算采用了基于ENDF／B－Ⅵ評價核數(shù)據(jù)庫的多群截面庫MATXS10，包括30群中子和12群光子的截面。

NUREG／CR－6115基準報告中，給出了采用MCNP／4A及DORT計算得到的軸向峰值處和焊縫處壓力容器內不同厚度處（內表面、1／4T）的快中子注量率（E＞1.0MeV）圓周方向分布結果。本文采用 MC－SN耦合程序計算，由于MATXS10庫中子能群在1.0MeV處并無能群邊界，因此對0.823～1.35MeV能群段內的中子進行插值處理，得到E＞1.0MeV的快中子注量率，并與基準報告提供的結果進行比較。分別如圖3～圖5所示。

圖2 b 1／8壓水堆模型垂直剖面圖Fig.2b Vertical section of 1／8PWR model

圖3 軸向峰值處壓力容器內表面快中子注量率（E＞1.0MeV）圓周分布曲線Fig.3 Circular distribution of E＞1.0MeV flux at pressure vessel inner wall axial peak location

圖3給出了軸向峰值處壓力容器內表面E＞1.0MeV的快中子注量率周向分布。從中可以看出，MC－SN計算得到的快中子注量率結果與NUREG／CR－6115基準報告中的 MCNP和DORT結果趨勢基本一致，吻合較好。與MCNP結果相比最大誤差小于7.45%，平均誤差小于0.25%。

圖4 軸向峰值處壓力容器1／4壁厚處快中子注量率（E＞1.0MeV）圓周分布曲線Fig.4 Circular distribution of E＞1.0MeV flux at pressure vessel T／4axial peak location

圖4給出了軸向峰值處壓力容器1／4壁厚處E＞1.0MeV的快中子注量率周向分布。從中可以看出，MC－SN計算得到的壓力容器快中子注量率結果與NUREG／CR－6115報告結果趨勢一致，與MCNP結果相比最大誤差小于11.96%，平均誤差為3.68%。

圖5 焊縫處壓力容器內表面快中子注量率（E＞1.0MeV）圓周分布Fig.5 Circular distribution of E＞1.0MeV flux at pressure vessel lower weld

圖5給出了焊縫處壓力容器內表面E＞1.0MeV的快中子注量率周向分布。從中可以看出，MC－SN計算得到的壓力容器快中子注量率結果與 NUREG／CR－6115PWR報告結果差別較小，與MCNP結果相比最大誤差小于7.57%，平均誤差為2.64%且趨勢一致。

從圖3～圖5可以看出，MC－SN計算結果與基準報告中MCNP和DORT計算結果相比基本吻合，少數(shù)角度對應的誤差相對大些，產(chǎn)生的原因主要有：（1）不同的計算程序基于的理論方法不同；（2）程序使用的核截面庫不同，MCNP程序使用的是連續(xù)截面數(shù)據(jù)庫，而TORT和DORT程序使用的是多群截面數(shù)據(jù)庫，理論精度低于連續(xù)截面數(shù)據(jù)庫；（3）幾何建模差別。MCNP較DORT程序對幾何模型和源項進行更精細的描述，減少了模型簡化近似引入的偏差，使其結果更可信。

3 結論

結合MC精確模擬復雜幾何和SN適合解決深穿透問題優(yōu)勢的三維MC－SN耦合方法，用于解決大型復雜核裝置屏蔽計算問題的難題。為了驗證耦合方法及程序系統(tǒng)的可靠性，采用BNL壓水堆基準例題對其進行了基準測試。耦合計算取得了與報告提供的與基準報告提供的MCNP和DORT基本一致的結果，驗證了方法的有效性和程序使用的正確性。表明三維MC－SN耦合方法可以為屏蔽系統(tǒng)優(yōu)化設計提供有力的技術支持。

[1]Zhang Bin，Chen Yixue，Wang Weijin，et al.Preliminary Shielding Analysis in Support of the CSNS Target Station Shutter Neutron Beam Stop Design[J］.CPC（HEP ＆NP），2011，35（8）：791－795.

[2]靳忠敏，陳義學，石生春，等 .基于MCNP的壓力容器快中子注量率計算參數(shù)敏感性分析[J］.原子能科學技術，2011，45（2）：195－199.

[3]李碩，陳義學，楊壽海 .基于二維離散縱標法的壓力容器快中子注量計算參數(shù)不確定度分析[C］.第十三屆反應堆數(shù)值計算與粒子輸運學術會議 .西安，2010.

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

[5]陳義學，吳宜燦，U.Fischer.國際聚變材料輻照裝置屏蔽中子學設計研究[J］.原子核物理評論，2006，23（2）：170－173.

[6]鐘兆鵬，施工，胡永明 .用 MCNP程序計算水平輻照孔道屏蔽 [J］.清華大學學報（自然科學版），2011，41（12）：16－18.

[7]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.

[8]Gabriel T A et al.，CALOR：A Monte Carlo Program Package for the Design and Analysis of Calorimeter Systems，ORNL／TM－5619[R］.USA：ORNL，1977.

[9]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.

[10]Urban W T，Seed T J，Dudziak D J.Nucleonic Analysis of a Preliminary Design for the ETF Neutral－beam－injector Duct Shielding[R］.LA－UR－80－2926.New Mexico：Los Alamos Scientific Laboratory，1980.

[11]Chen Y.Coupled Monte Carlo－discrete Ordinates Computational Scheme for Three－dimensional Shielding Calculations of Large and Complex Nuclear Facilities[D］.Germany：Forschungszentrum Karlsruhe，2005.

[12]Y.Chen， U.Fischer.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，2005，A（551）：387－395.

[13]BRIESMEISTER J F.MCNP：A general Monte Carlo N－particle Transport Code，Version 4C，LA－13709－M[R］.USA：Los Alamos National Laboratory，2000.

[14]RhoadesW A，Simpson D B.The TORT Three－Dimensional Discrete Ordinates Neutron／Photon Trans Port Code[R］.ORNL／TM13221，Oak Ridge National Laboratory，1997.

[15]Carew J F，Hu K，Aroson A，et al.PWR and BWR Pressure Vessel Fluence Calculation Benchmark Problem and Solutions[R］.Washington：Brookhaven National Laboratory，2001.