韓靜茹，陳義學(xué)，袁龍軍，張春明

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

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

1 三維SN-MC耦合系統(tǒng)框架

TDOMINO耦合系統(tǒng)流程圖示于圖1。以現(xiàn)有程序?yàn)榛A(chǔ)，通過(guò)開發(fā)建立相應(yīng)的耦合接口和源粒子抽樣程序來(lái)實(shí)現(xiàn)程序數(shù)據(jù)交換和計(jì)算。這樣的耦合框架既避免了對(duì)現(xiàn)有程序結(jié)構(gòu)的修改，也便于今后對(duì)程序進(jìn)行維護(hù)和升級(jí)。同時(shí)也使耦合形式具有較好的靈活性，可根據(jù)問(wèn)題分析的需要選擇用于耦合計(jì)算的程序。構(gòu)成TDOMINO耦合系統(tǒng)的基礎(chǔ)程序包括：三維離散縱標(biāo)程序TORT[5]和蒙特卡羅程序MCNP[6]。在耦合系統(tǒng)中，為充分發(fā)揮兩種程序的優(yōu)勢(shì)，TORT和MCNP分別用以模擬大塊屏蔽區(qū)和復(fù)雜幾何區(qū)計(jì)算。其中三維SN計(jì)算得到連接面的粒子角通量密度，并以BNDRYS格式存儲(chǔ)輸出，作為SN-MC接口程序的輸入文件之一。SN-MC接口程序?qū)⑦B接面的角通量密度轉(zhuǎn)換為粒子概率分布，再經(jīng)MC-SOURCE源粒子抽樣程序轉(zhuǎn)換為MC源粒子信息參數(shù)，為MC計(jì)算提供源項(xiàng)，實(shí)現(xiàn)三維SN和MC耦合輸運(yùn)計(jì)算，得到復(fù)雜幾何區(qū)的粒子通量分布。該程序系統(tǒng)可處理三維x-y-z及r-θ-z幾何。

圖1 TDOMINO耦合系統(tǒng)流程圖

2 三維SN-MC耦合方法

對(duì)于上述兩個(gè)基礎(chǔ)程序，由于其物理建模和數(shù)值求解方法不同，耦合需解決連接區(qū)內(nèi)SN計(jì)算的粒子角通量密度分布和MC計(jì)算所需源項(xiàng)之間映射的問(wèn)題。在TDOMINO耦合系統(tǒng)中，通過(guò)開發(fā)三維SN-MC接口程序和MC-SOURCE源粒子抽樣程序分別實(shí)現(xiàn)連接區(qū)內(nèi)SN計(jì)算得到的粒子角通量密度和粒子概率分布之間映射及粒子概率分布和MC源粒子信息間的轉(zhuǎn)換。具體算法為：

(1)

(2)

(3)

(4)

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

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

3 耦合程序的驗(yàn)證

為驗(yàn)證三維SN-MC耦合程序，采用美國(guó)橡樹嶺國(guó)家實(shí)驗(yàn)室公布的HBR-2屏蔽基準(zhǔn)題[7]，利用TDOMINO耦合系統(tǒng)進(jìn)行驗(yàn)證計(jì)算，并針對(duì)圓柱坐標(biāo)系和直角坐標(biāo)系下的耦合模型特點(diǎn)，采用不同的耦合形式進(jìn)行計(jì)算。

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

3.1 HBR-2基準(zhǔn)題

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

3.2 計(jì)算模型

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

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

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

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

a——耦合計(jì)算模型；b——TORT計(jì)算模型；c——MCNP計(jì)算模型

3.3 計(jì)算結(jié)果

分別采用三維SN-MC耦合程序TDOMINO和單獨(dú)的MCNP程序?qū)BR-2基準(zhǔn)模型輻照監(jiān)督管處6種典型核素的比活度進(jìn)行計(jì)算，并與實(shí)驗(yàn)測(cè)量、離散縱標(biāo)法程序計(jì)算結(jié)果進(jìn)行對(duì)比分析。

表1列出HBR-2基準(zhǔn)模型輻照監(jiān)督管處6種典型核素比活度的實(shí)驗(yàn)測(cè)量值。表2列出不同程序計(jì)算的6種核素比活度與實(shí)驗(yàn)測(cè)量值間的比值(C/M)。其中實(shí)驗(yàn)測(cè)量值和DORT采用BUGLE-96[8]庫(kù)的計(jì)算結(jié)果取自HBR-2基準(zhǔn)報(bào)告，TORT結(jié)合TEXT10多群截面數(shù)據(jù)庫(kù)的計(jì)算結(jié)果取自文獻(xiàn)[9]，SN-MC和SN-MC-2分別表示TDOMINO在圓柱坐標(biāo)系和直角坐標(biāo)系下的計(jì)算結(jié)果。MCNP表示堆芯組件均勻模型的單一MC計(jì)算結(jié)果，MCNP-PIN表示堆芯PIN功率的計(jì)算結(jié)果取自文獻(xiàn)[10]。DORT、SN-MC、SN-MC-2、MCNP、MCNP-PIN和TORT程序計(jì)算得到的輻照監(jiān)督管處的平均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。計(jì)算結(jié)果表明，TDOMINO耦合系統(tǒng)計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量值和其他程序計(jì)算值吻合較好，初步驗(yàn)證了三維SN-MC耦合系統(tǒng)TDOMINO在不同空間坐標(biāo)系中應(yīng)用的正確性。在上述計(jì)算中，SN-MC耦合程序的計(jì)算時(shí)間為3 908 min，其中，SN-MC中MCNP計(jì)算耗時(shí)3 869 min，最大統(tǒng)計(jì)誤差為3%左右，SN-MC中TORT計(jì)算耗時(shí)39 min。單一MCNP程序計(jì)算耗時(shí)7 180 min，最大統(tǒng)計(jì)誤差為5%左右。表明三維SN-MC耦合程序與單一MCNP程序相比，在較少的時(shí)間下得到的統(tǒng)計(jì)誤差更小，計(jì)算結(jié)果可靠性更強(qiáng)。

表1 實(shí)驗(yàn)測(cè)量所得的比活度

表2 計(jì)算與測(cè)量的比活度之比

4 結(jié)論

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

參考文獻(xiàn)：

[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] 鄭征，吳宏春，曹良志，等. 蒙特卡羅與離散縱標(biāo)耦合方法在壓水堆堆腔漏束計(jì)算中的應(yīng)用[J]. 強(qiáng)激光與粒子束,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] 楊壽海,陳義學(xué),王偉金,等. 三維離散縱標(biāo)方法在RPV快中子注量率計(jì)算中的初步應(yīng)用[J]. 核科學(xué)與工程，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] 石生春. 基于蒙特卡羅方法的壓水堆壓力容器快中子注量率的計(jì)算分析[D]. 北京：華北電力大學(xué)，2010.