韓靜茹，陳義學(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ì)算方法是保證屏蔽系統(tǒng)設(shè)計(jì)質(zhì)量的關(guān)鍵。華北電力大學(xué)在輻射屏蔽方法及應(yīng)用方面開(kāi)展了大量的研究工作[1－5］。蒙特卡羅方法（MC）和離散縱標(biāo)方法（SN）是最常用的屏蔽計(jì)算方法。MC方法的優(yōu)點(diǎn)是可以精確模擬復(fù)雜幾何模型，但計(jì)算耗時(shí)，尤其對(duì)于屏蔽計(jì)算中常見(jiàn)的深穿透問(wèn)題，有的甚至根本無(wú)法得出計(jì)算結(jié)果。盡管有MC分段計(jì)算[6］方法的研究用于解決深穿透的問(wèn)題，但由于分段計(jì)算中MC統(tǒng)計(jì)誤差的連續(xù)傳遞性，但這仍舊是其在研究和實(shí)際應(yīng)用中最大的困難。與MC方法相比，SN方法特別適合解決深穿透問(wèn)題但卻難以精確描述復(fù)雜幾何模型。對(duì)于大型復(fù)雜核裝置屏蔽計(jì)算問(wèn)題，MC方法或SN無(wú)法提供可靠的計(jì)算結(jié)果。

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

1 MC－SN耦合方法及程序

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

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

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

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

2 MC－SN基準(zhǔn)測(cè)試

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

2.1 基準(zhǔn)模型描述

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

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

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

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

2.2 計(jì)算結(jié)果及分析

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

[5]陳義學(xué)，吳宜燦，U.Fischer.國(guó)際聚變材料輻照裝置屏蔽中子學(xué)設(shè)計(jì)研究[J］.原子核物理評(píng)論，2006，23（2）：170－173.

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