基于蒙特卡羅-離散縱標(biāo)雙向耦合方法的堆坑粒子注量率計(jì)算

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

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

在核電站設(shè)計(jì)中，反應(yīng)堆壓力容器位于反應(yīng)堆堆坑內(nèi)，壓力容器和一次混凝土屏蔽墻之間構(gòu)成的空間一般稱為堆腔。反應(yīng)堆運(yùn)行過(guò)程中，堆芯中的一部分中子和光子從壓力容器外表面逸出，進(jìn)入堆腔，并繼續(xù)向外泄漏，有一部分泄漏到壓力容器下方的堆坑底部。堆坑底部設(shè)有堆坑通道、小室及堆內(nèi)測(cè)量?jī)x表間，用于壓力容器絕熱層、壓力容器監(jiān)測(cè)設(shè)備及堆內(nèi)測(cè)量?jī)x表的安裝和維修。堆坑粒子注量率計(jì)算分析對(duì)于核電站設(shè)計(jì)和安全分析具有重要意義。

堆坑底部粒子注量率計(jì)算涉及的橫向、縱向尺寸均達(dá)十多米，且中間的屏蔽構(gòu)件幾何復(fù)雜，顯然這是一大尺寸復(fù)雜幾何系統(tǒng)內(nèi)的輸運(yùn)問(wèn)題。由于反應(yīng)堆廠房幾何尺寸較大，且堆坑通道軸線與反應(yīng)堆主軸垂直，則自壓力容器外表面泄漏的粒子到達(dá)小室屏蔽門附近的幾率很小。直接采用單一的蒙特卡羅(MC)方法難以在有效的時(shí)間內(nèi)解決此類小概率事件問(wèn)題。而在壓力容器外，由于主冷卻劑進(jìn)、出口接管，堆外探測(cè)器孔道，堆坑通道及小室等的存在，堆腔幾何結(jié)構(gòu)極為復(fù)雜，單一的離散縱標(biāo)(SN)法難以精確描述堆坑及小室的復(fù)雜幾何，使計(jì)算結(jié)果產(chǎn)生較大誤差。因此，本文提出結(jié)合兩種方法優(yōu)勢(shì)的三維蒙特卡羅-離散縱標(biāo)雙向耦合(MC-SN雙向耦合)方法，用于壓水堆堆坑底部輻射泄漏粒子注量率的計(jì)算。

1 MC-SN雙向耦合方法

MC和SN方法是求解中子輸運(yùn)問(wèn)題的兩種不同方法。MC方法通過(guò)對(duì)大量單個(gè)粒子的運(yùn)動(dòng)歷史逐個(gè)進(jìn)行跟蹤模擬，然后用統(tǒng)計(jì)方法給出隨機(jī)變量某個(gè)特征的估計(jì)量，作為問(wèn)題的解，SN方法是在相空間(r×E×Ω)內(nèi)對(duì)空間r、能量E和方向變量Ω采用直接離散方法數(shù)值求解中子輸運(yùn)方程，得到網(wǎng)格ri、能群Eg和離散方向Ωm區(qū)間內(nèi)的粒子角通量密度ψ(ri,Eg,Ωm)。實(shí)現(xiàn)蒙特卡羅方法和離散縱標(biāo)方法之間的耦合計(jì)算，關(guān)鍵在于MC模擬的粒子徑跡信息與SN方法計(jì)算的角通量密度之間的相互轉(zhuǎn)換。其中，MC粒子徑跡信息到SN角通量密度的轉(zhuǎn)換關(guān)系[1]如下：

(1)

式中：ψ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)系[2]如下：

(2)

(3)

(4)

式中：λm為粒子飛行方向與面法線方向夾角的余弦值；Ai為網(wǎng)格區(qū)間ri對(duì)應(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)的概率。

2 MC-SN雙向耦合系統(tǒng)框架

基于蒙特卡羅-離散縱標(biāo)法的雙向耦合計(jì)算框架如圖1所示。根據(jù)MC-SN雙向耦合轉(zhuǎn)換方法，通過(guò)自主研發(fā)接口程序來(lái)實(shí)現(xiàn)。三維MC-SN雙向耦合系統(tǒng)流程如圖2所示，其中MC程序采用國(guó)際通用的粒子輸運(yùn)程序MCNP[3]；SN程序采用美國(guó)橡樹嶺國(guó)家實(shí)驗(yàn)室開發(fā)的三維離散縱標(biāo)程序TORT[4]。在耦合系統(tǒng)中，為了充分發(fā)揮兩種方法的優(yōu)勢(shì)，MCNP和TORT分別用以模擬復(fù)雜幾何區(qū)和大塊屏蔽區(qū)的計(jì)算。其中MCNP計(jì)算得到的穿過(guò)指定面的粒子數(shù)據(jù)被精確記錄在SSW面源文件中，通過(guò)MC-SN接口程序用于生成TORT程序的邊界源文件，作為TORT程序的輸入之一，實(shí)現(xiàn)三維MC和SN耦合輸運(yùn)計(jì)算；TORT計(jì)算得到連接面的粒子角通量密度，SN-MC接口程序?qū)⑦B接面的角通量密度轉(zhuǎn)換為粒子概率分布，再經(jīng)MC-SOURCE源粒子抽樣程序轉(zhuǎn)換為MC源粒子信息參數(shù)，為MCNP計(jì)算提供源項(xiàng)，實(shí)現(xiàn)三維SN和MC耦合輸運(yùn)計(jì)算。

圖1 MC-SN雙向耦合系統(tǒng)框架

MC-SN雙向耦合系統(tǒng)可處理三維x-y-z及R-θ-z幾何。對(duì)于大尺寸復(fù)雜幾何系統(tǒng)屏蔽計(jì)算問(wèn)題，可選取MC-SN雙向耦合方法，針對(duì)裝置自身特點(diǎn)，設(shè)計(jì)較優(yōu)的分割方案，靈活選擇不同的耦合方法和程序進(jìn)行屏蔽計(jì)算。既可提高對(duì)大型復(fù)雜核裝置屏蔽問(wèn)題計(jì)算結(jié)果的可靠性和精確度，同時(shí)又提高了收斂速度，節(jié)省了機(jī)時(shí)。

3 計(jì)算模型及步驟

3.1 計(jì)算模型

某反應(yīng)堆及堆坑布置如圖3所示。堆坑通道和小室相連通，布置在51°方位上，與反應(yīng)堆主軸線相垂直。屏蔽門在小室側(cè)部，通風(fēng)孔布置在小室頂部[5]。根據(jù)計(jì)算模型特點(diǎn)，采用MC-SN-MC耦合程序系統(tǒng)對(duì)堆坑底部通道和小室進(jìn)行中子、光子注量率分布計(jì)算，用以提供小室屏蔽門設(shè)計(jì)所需的輻照強(qiáng)度數(shù)據(jù)。

圖2 三維MC-SN雙向耦合系統(tǒng)流程圖

a——垂直剖面圖；b——水平剖面圖；c——堆坑底部及探測(cè)器示意圖

本計(jì)算模型中，系統(tǒng)坐標(biāo)原點(diǎn)設(shè)置在堆芯中平面處，坐標(biāo)系z(mì)軸與反應(yīng)堆軸心線重合，x軸和y軸分別與0°方位和90°方位重合。計(jì)算中設(shè)置了兩個(gè)耦合面S1和S2(圖3a)。位于堆芯下表面上的S1，其為高度z=-182.88 cm、半徑為221.655 cm的90°圓平面，在此表面上實(shí)現(xiàn)MC-SN耦合計(jì)算；位于壓力容器底面，高度z=-413.526 cm處，半徑為221.655 cm的90°圓平面S2為第2耦合面，在此表面上實(shí)現(xiàn)SN-MC耦合。P1、P2和P3分別表示反應(yīng)堆堆坑通道、小室及小室屏蔽門附近的探測(cè)器(圖3c)，均為半徑為18 cm的小球。表1列出3個(gè)探測(cè)器的具體坐標(biāo)位置。

表1 反應(yīng)堆堆坑探測(cè)器的位置

3.2 計(jì)算步驟

三維MC-SN-MC耦合計(jì)算具體步驟如下。1) MC計(jì)算：采用MCNP建立反應(yīng)堆下反射層以上區(qū)域的精細(xì)模型進(jìn)行計(jì)算，求得穿過(guò)S1表面的粒子(中子、光子)徑跡信息，并存儲(chǔ)到WSSA文件。2) MC-SN計(jì)算：采用MC-SN接口程序?qū)SSA文件中記錄的粒子徑跡信息轉(zhuǎn)換為角通量密度分布，并生成SN邊界源文件。3) SN計(jì)算：采用TORT建立S1表面到壓力容器底面的R-θ-z模型，選取P3階勒讓德展開，S8全對(duì)稱高斯求積組，進(jìn)行邊界源計(jì)算，得到S2面粒子角通量密度分布，完成三維MC-SN耦合計(jì)算。4) SN-MC計(jì)算：通過(guò)SN-MC接口程序，將TORT計(jì)算得到的S2面上的粒子角通量密度分布轉(zhuǎn)換為中子在空間、能量和方向上的概率分布函數(shù)。5) MC計(jì)算：在壓力容器底面S2以下，采用MCNP建立堆坑精細(xì)模型，通過(guò)已嵌入MCNP程序的源抽樣程序MC-SOURCE對(duì)第4步中的概率分布函數(shù)進(jìn)行抽樣，跟蹤模擬S2面上的粒子，計(jì)算求得堆坑底部3個(gè)探測(cè)器的中子和光子注量率分布，完成三維SN-MC耦合計(jì)算。耦合計(jì)算中MCNP計(jì)算采用點(diǎn)截面數(shù)據(jù)庫(kù)，TORT計(jì)算采用包含30群中子和12群光子的TEXT10多群數(shù)據(jù)庫(kù)。本計(jì)算中的模型結(jié)構(gòu)尺寸及材料與實(shí)際反應(yīng)堆存在一定差別，計(jì)算結(jié)果為1個(gè)中子源的歸一化結(jié)果。

4 計(jì)算結(jié)果分析及討論

4.1 MC-SN-MC耦合計(jì)算結(jié)果

采用三維MCNP-TORT-MCNP耦合程序，設(shè)置F4體探測(cè)器計(jì)數(shù)，計(jì)算得到3個(gè)探測(cè)器處的中子和光子注量率。同時(shí)根據(jù)NCRP-38[6]中的中子通量-劑量轉(zhuǎn)換因子和光子通量-劑量轉(zhuǎn)換因子求得指定探測(cè)器處的劑量率。耦合計(jì)算中，最耗時(shí)的為上節(jié)計(jì)算步驟中所描述的第1步和第5步中的MCNP程序計(jì)算過(guò)程。其中第1步MCNP計(jì)算300 min，S1面獲得1.0×106個(gè)粒子徑跡，第5步中MCNP中子光子輸運(yùn)計(jì)算分別進(jìn)行，其中中子計(jì)算400 min，光子計(jì)算2 000 min，獲得P1、P2和P3處中子和光子注量率分布。MC-SN-MC計(jì)算得到的3個(gè)探測(cè)器處歸一化的中子注量率、中子劑量率、光子注量率及光子劑量率值列于表2?？煽闯觯憾芽拥撞客ǖ馈⑿∈壹靶∈移帘伍T附近，對(duì)輻射劑量安全起主導(dǎo)作用的是中子，光子對(duì)劑量率的貢獻(xiàn)較中子小1～2個(gè)數(shù)量級(jí)。

3個(gè)探測(cè)器處的中子能譜和光子能譜示于圖4?？煽闯?，P1、P2和P3處的粒子能譜分布趨勢(shì)基本一致，小室屏蔽門附近P3探測(cè)器處的中子和光子注量率均為最低。

表2 堆坑底部探測(cè)器處粒子注量率和劑量率

注：括號(hào)內(nèi)為統(tǒng)計(jì)誤差

圖4 探測(cè)器處中子能譜和光子能譜

4.2 計(jì)算結(jié)果驗(yàn)證比較

采用單一的MCNP程序建立整個(gè)模型，并采取相同的探測(cè)器計(jì)數(shù)，計(jì)算34 000 min后，跟蹤模擬2.0×108個(gè)粒子，P1、P2和P3探測(cè)器處中子注量率的統(tǒng)計(jì)誤差分別達(dá)到4.52%、7.78%和10.37%。由此可見(jiàn)，單一的蒙特卡羅方法在求解深穿透復(fù)雜屏蔽計(jì)算問(wèn)題時(shí)很難收斂。選取單一的MCNP計(jì)算相對(duì)統(tǒng)計(jì)誤差最小的P1探測(cè)器處的中子注量率分布與三維MCNP-TORT-MCNP計(jì)算結(jié)果比較，結(jié)果如圖5所示，中子能譜分布趨勢(shì)相同。耦合程序計(jì)算的P1處的中子注量率較MCNP的計(jì)算結(jié)果偏低約2.92%，主要是由于MCNP程序計(jì)算結(jié)果相比耦合計(jì)算統(tǒng)計(jì)誤差相對(duì)偏大。單一MCNP程序隨著計(jì)算時(shí)間的增加，結(jié)果呈變小趨勢(shì)。初步證明了三維MC-SN雙向耦合程序系統(tǒng)的正確性和有效性。

圖5 P1處中子能譜

5 小結(jié)

為解決大尺寸復(fù)雜幾何系統(tǒng)輸運(yùn)問(wèn)題難以采用單一的模型及程序進(jìn)行求解的難題，開發(fā)了三維蒙特卡羅-離散縱標(biāo)雙向耦合方法。將基于此耦合方法的MCNP-TORT-MCNP程序用于某壓水堆堆坑粒子注量率計(jì)算，并與單一的MCNP程序計(jì)算進(jìn)行比較，結(jié)果表明，蒙特卡羅-離散縱標(biāo)雙向耦合程序系統(tǒng)在速度和精度上均有很大優(yōu)勢(shì)，驗(yàn)證了三維蒙特卡羅-離散縱標(biāo)雙向耦合方法的正確性和可靠性。采用蒙特卡羅-離散縱標(biāo)雙向耦合程序系統(tǒng)使大型復(fù)雜核裝置全三維屏蔽計(jì)算成為可能。

參考文獻(xiàn)：

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

[2]韓靜茹，陳義學(xué)，石生春，等. 基于離散縱標(biāo)法與蒙特卡羅方法的三維耦合程序開發(fā)[J]. 核科學(xué)與工程，2012，32(2)：160-164.

HAN Jingru, CHEN Yixue, SHI Shengchun, et al. Development of three dimensional coupled code system based on discrete ordinates and Monte Carlo method[J]. Nuclear Science and Engineering, 2012, 32(2): 160-164(in Chinese).

[3]BRIESMEISTER J F. MCNP: A general Monte CarloN-particle transport code, version 4C, LA-13709-M[R]. USA: Los Alamos National Laboratory, 2000.

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

[5]劉桂蓮. 秦山核電二期工程反應(yīng)堆堆坑底部輻射通量分布計(jì)算[J]. 核動(dòng)力工程，1998，19(4)：375-379.

LIU Guilian. Radiation flux calculation of reactor pit bottom in Qinshan Phase Ⅱ Nuclear Power project[J]. Nuclear Power Engineering, 1998, 19(4): 375-379(in Chinese).

[6]ROSSI H H. Protection against neutron radiation, NCRP-38[R]. US: National Council on Radiation Protection and Measurement, 1971.