節(jié)塊法與蒙特卡羅方法耦合計(jì)算研究

王新哲，喻 宏，徐 李，胡 赟，楊曉燕

（中國(guó)原子能科學(xué)研究院快堆研究設(shè)計(jì)所，北京 102413）

節(jié)塊法與蒙特卡羅方法耦合計(jì)算研究

王新哲，喻 宏，徐 李，胡 赟，楊曉燕

（中國(guó)原子能科學(xué)研究院快堆研究設(shè)計(jì)所，北京 102413）

普通節(jié)塊法無(wú)法在計(jì)算中獲得不同組件內(nèi)精細(xì)中子通量密度分布的信息。本文提出一種利用入射角通量將節(jié)塊法與蒙特卡羅方法相耦合的方法（節(jié)塊－蒙卡入射角通量耦合方法），并編制了計(jì)算程序進(jìn)行驗(yàn)證。結(jié)果表明：本文計(jì)算結(jié)果與參考值相符，節(jié)塊－蒙卡入射角通量耦合方法適用于局部特定位置精細(xì)中子通量密度等參數(shù)的計(jì)算，計(jì)算效率高，計(jì)算結(jié)果準(zhǔn)確。

節(jié)塊法；蒙特卡羅方法；耦合；入射角通量；精細(xì)中子通量密度

堆芯精細(xì)中子通量密度和功率分布的計(jì)算是反應(yīng)堆中子學(xué)分析軟件發(fā)展的方向之一，準(zhǔn)確獲得堆芯精細(xì)中子通量密度和功率分布，可有效減少設(shè)計(jì)中的保守量，提高反應(yīng)堆的經(jīng)濟(jì)性。普通節(jié)塊法在計(jì)算中無(wú)法獲得不同組件內(nèi)精細(xì)中子通量密度分布等參數(shù)的信息?，F(xiàn)階段，精細(xì)中子通量密度與功率分布計(jì)算的主要方法有功率重構(gòu)方法［1］、蒙特卡羅方法［2］及特征線方法［3］等。這些方法在現(xiàn)有條件下皆存在一定局限性，如功率重構(gòu)方法對(duì)于復(fù)雜結(jié)構(gòu)組件計(jì)算精度差［4］，蒙特卡羅方法與特征線方法的計(jì)算速度較慢。

本文通過(guò)分析節(jié)塊法與蒙特卡羅方法的特點(diǎn)，提出一種利用入射角通量將二者相耦合的方法——節(jié)塊－蒙卡入射角通量耦合方法，對(duì)任意幾何組件的精細(xì)中子通量密度與功率分布進(jìn)行計(jì)算。

1 理論模型

1.1 節(jié)塊－蒙卡入射角通量耦合方法描述

節(jié)塊－蒙卡入射角通量耦合方法的基本思想是，將求解中子輸運(yùn)方程的兩種思路——確定論方法與蒙特卡羅方法通過(guò)邊界流數(shù)據(jù)傳遞的方式耦合在一起，從而發(fā)揮二者各自的優(yōu)勢(shì)。

穩(wěn)態(tài)中子輸運(yùn)方程為：

其中：Φ為中子通量密度；r為中子空間位置向量；Ω、Ω′分別為碰撞前后中子運(yùn)動(dòng)方向的單位向量；E、E′分別為碰撞前后中子的能量；Σt為材料總宏觀截面；Σs為材料散射截面；f為散射函數(shù)；Qf為裂變?cè)错?xiàng)；S為外源項(xiàng)。

確定論方法的基本思想是離散，中子輸運(yùn)方程中的未知參數(shù)中子通量密度Φ（r，E，Ω）是空間位置、能量與飛行角度的函數(shù)，這些自變量在物理上都是連續(xù)的。為了將連續(xù)的中子輸運(yùn)方程轉(zhuǎn)換為有限個(gè)可解的數(shù)學(xué)物理方程，必須對(duì)這些連續(xù)參數(shù)進(jìn)行離散。現(xiàn)有確定論方法對(duì)于幾何空間的離散，大多采用規(guī)則網(wǎng)格劃分，不能進(jìn)行精細(xì)幾何描述，對(duì)于結(jié)構(gòu)復(fù)雜的組件設(shè)計(jì)有一定計(jì)算誤差。

蒙特卡羅方法能進(jìn)行精細(xì)幾何描述，問(wèn)題適應(yīng)能力強(qiáng)，但由于其需對(duì)大量粒子進(jìn)行統(tǒng)計(jì)，因此全堆計(jì)算需較長(zhǎng)時(shí)間。

在某些計(jì)算中只關(guān)注堆內(nèi)特定位置處的精細(xì)分布，上述兩種方法很難快速、準(zhǔn)確地給出計(jì)算結(jié)果。此時(shí)，可利用將二者耦合的思路，具體計(jì)算流程如下：1）利用節(jié)塊法程序進(jìn)行快速全堆計(jì)算，得出中子通量密度分布和有效增殖因數(shù)等信息；2）將獲得的信息轉(zhuǎn)換成所求關(guān)注位置處組件的入射角通量；3）利用獲得的入射角通量信息產(chǎn)生蒙特卡羅程序的面源文件，同時(shí)按照精細(xì)幾何進(jìn)行組件描述，進(jìn)行蒙特卡羅面源計(jì)算。

通過(guò)上述3步，即可以較低的計(jì)算代價(jià)獲得特定組件的精細(xì)幾何描述計(jì)算參數(shù)。

1.2 入射角通量分布計(jì)算

本文采用六角形節(jié)塊程序NAS［5］作為節(jié)塊法程序。NAS由中國(guó)原子能科學(xué)研究院快堆研究設(shè)計(jì)所開(kāi)發(fā)，能按照擴(kuò)散與輸運(yùn)兩種方法進(jìn)行三維堆芯計(jì)算，具有微擾計(jì)算、時(shí)空動(dòng)力學(xué)計(jì)算、燃耗計(jì)算、堆芯燃耗管理及優(yōu)化等功能。NAS在徑向上按照組件劃分節(jié)塊，軸向則可根據(jù)堆芯布置適當(dāng)進(jìn)行劃分。本文堆芯計(jì)算中，采用NAS輸運(yùn)模塊進(jìn)行，表面通量采用雙球諧DP3近似展開(kāi)，節(jié)塊的求解采用了橫向積分方法，它將三維輸運(yùn)方程轉(zhuǎn)化為4個(gè)一維輸運(yùn)方程，并采用一維離散縱標(biāo)（SN）方法進(jìn)行求解。

NAS通過(guò)堆芯計(jì)算，可給出節(jié)塊各面上的入射通量Φin、入射流Jin、二階入射通量Φin2及三階入射通量Φin3，節(jié)塊間的耦合關(guān)系為：

其中：j′和j分別為臨近的兩個(gè)節(jié)塊；ψij為j節(jié)塊i方向的入射角通量；μi和ωi分別為SN方法中高斯－勒讓德求積組的離散方向余弦與求積權(quán)重，其中μi轉(zhuǎn)換到［0，1］區(qū)間，∑iωi＝1；Pn為n階勒讓德多項(xiàng)式。在現(xiàn)版本NAS中，角度的離散只考慮了與節(jié)塊表面法線夾角一個(gè)方向。

據(jù)此可得到入射角通量的表達(dá)式：

其中：

通過(guò)式（3）與式（4）即可求得各面上分能群入射中子角通量ψ（Si，Ej，μl），其中，Si為第i面，Ej為第j能群，μl為第l離散角度。

要獲得蒙特卡羅程序所需的面源分布，需將ψ（Si，Ej，μl）轉(zhuǎn)換為面入射粒子數(shù)，按照式（5）進(jìn)行轉(zhuǎn)換。

其中：N（Si，Ej，μl）為面入射粒子數(shù)；Ai為入射面面積。

對(duì)于總計(jì)算粒子數(shù)NT，各面上、各能群、各角度的入射中子數(shù)N′（Si，Ej，μl）有：

1.3 蒙特卡羅面源抽樣方法

本文中蒙特卡羅程序采用MCNP［6］，MCNP可讀取面源文件RSSA進(jìn)行輸運(yùn)計(jì)算。本文根據(jù)式（6）中的中子數(shù)分布通過(guò)蒙特卡羅方法進(jìn)行抽樣，自動(dòng)生成該文件。面源文件需提供每個(gè)中子的編號(hào)、時(shí)間、權(quán)重、空間坐標(biāo)、飛行方向、與入射面夾角和能量等信息，其中需要抽樣的信息為空間坐標(biāo)、飛行方向和能量。

由于NAS尚不能提供入射中子面分布信息，故在抽樣過(guò)程中認(rèn)為每個(gè)入射面上中子的分布為均勻分布。為方便抽樣，抽樣方法為先在一個(gè)標(biāo)準(zhǔn)矩形面上進(jìn)行抽樣，然后再將粒子坐標(biāo)轉(zhuǎn)換到相應(yīng)的面上。

對(duì)于角度分布，認(rèn)為在每個(gè)離散角度內(nèi)，中子飛行方向均勻分布，即：

其中，f（φ（μ））為飛行方向?yàn)棣痰闹凶油棵芏鹊母怕拭芏取?/p>

則入射粒子與面夾角的概率密度函數(shù)為：

其中，C為歸一化常數(shù)。

由式（7）、（8）可見(jiàn)，入射粒子關(guān)于面夾角的分布為一次分布。抽樣過(guò)程中也需進(jìn)行相對(duì)面位置的轉(zhuǎn)換，在圖1所示編號(hào)為1的面上，按照各向同性空間粒子抽樣［7］，然后按式（8）以一定概率對(duì)其取舍，獲得夾角參數(shù)sinθ、cosθ、sinφ、cosφ，其中θ為極角，φ為方位角，則由式（9）可得到粒子的飛行方向。

圖1 六角形節(jié)塊面編號(hào)示意圖Fig.1 Scheme of surface number for hexagonal nodal

其中，U、V、W分別為中子與x、y、z軸夾角的余弦。

按照式（10），對(duì)上述夾角進(jìn)行旋轉(zhuǎn)編號(hào)，即可得到粒子飛行方向。

式中，N為面編號(hào)。

六角形節(jié)塊上頂面、下底面的入射中子分布可很容易地抽樣得到。

各能群內(nèi)能量的分布由程序輸入給出，可為裂變譜、1／E譜、麥克斯韋譜或均勻分布譜。其中，裂變譜采用乘減抽樣方法，1／E譜和麥克斯韋譜采用乘抽樣方法。本文中NAS采用171群的NVITAMIN－C庫(kù)，為方便快堆計(jì)算，將最后11群并為一群，每群內(nèi)能量分布按照NJOY程序手冊(cè)中用于VITAMIN－E庫(kù)制作的權(quán)重函數(shù)給出，并將其中高能部分的聚變譜改為1／E譜［8］。

2 程序開(kāi)發(fā)

按照上述方法，基于NAS與MCNP編制了用于六角形組件計(jì)算的MCNTC程序。MCNTC程序流程圖如圖2所示。

圖2 MCNTC程序流程圖Fig.2 Flow chart of MCNTC

由圖2可見(jiàn)，只需給出NAS輸入文件、MCNP組件描述文件及部分控制參數(shù)，即可完成程序功能。通過(guò)輸入可控制抽樣的能群劃分及群內(nèi)中子能譜、中子飛行方向離散角度劃分等參數(shù)。

3 程序驗(yàn)證

采用3個(gè)算例對(duì)節(jié)塊－蒙卡入射角通量耦合方法及MCNTC進(jìn)行驗(yàn)證。

3.1 中國(guó)實(shí)驗(yàn)快堆首爐堆芯燃料組件內(nèi)功率分布驗(yàn)證

中國(guó)實(shí)驗(yàn)快堆（CEFR）是我國(guó)第一座鈉冷快中子反應(yīng)堆，設(shè)計(jì)熱功率為65MW，于2010年7月首次達(dá)到臨界。在初裝料時(shí)，堆芯裝有79盒燃料組件、8個(gè)控制棒組件和1個(gè)中子源組件，堆芯外面為鋼反射層，共有336個(gè)反射層組件。反射層外為屏蔽層，裝載有230個(gè)硼屏蔽層組件［9］。燃料組件每盒有61根燃料棒，燃料芯塊為富集度64.4%的UO2。堆芯燃料組件編號(hào)如圖3所示。

圖3 堆芯燃料組件編號(hào)示意圖Fig.3 Scheme of fuel assembly number for reactor core

按照CEFR物理啟動(dòng)實(shí)驗(yàn)測(cè)得的臨界棒位，MCNTC中采用NAS輸運(yùn)功能進(jìn)行堆芯計(jì)算，keff的計(jì)算值為1.004 03，與實(shí)驗(yàn)值符合較好。耦合計(jì)算中MCNP計(jì)算采用燃料組件活性區(qū)精細(xì)描述，并進(jìn)行裂變功率統(tǒng)計(jì)，每次計(jì)算抽樣107個(gè)面源粒子，棒功率統(tǒng)計(jì)不確定度小于0.25%。同時(shí)，采用MCNP全堆計(jì)算作為參考結(jié)果，MCNP全堆計(jì)算中，燃料活性區(qū)為精細(xì)幾何描述，燃料其他區(qū)域及其他組件采用均勻幾何描述，并統(tǒng)計(jì)所有燃料棒的裂變功率，有效統(tǒng)計(jì)中子數(shù)為108，keff的計(jì)算值為0.999 91±0.000 07，棒功率統(tǒng)計(jì)不確定度小于0.20%。

MCNTC與MCNP皆在Intel Xeon E7－4850上運(yùn)行。MCNTC采用串行計(jì)算，每次程序運(yùn)行時(shí)間約為30min，其中約10min為堆芯臨界計(jì)算，其余約20min為面源生成及組件精細(xì)功率分布計(jì)算，79盒燃料組件全部計(jì)算約需40CPU小時(shí)。MCNP采用10核20線程并行計(jì)算，程序運(yùn)行時(shí)間約為2 500CPU小時(shí)。因此，MCNTC在全堆功率分布計(jì)算上的速度較傳統(tǒng)MCNP快數(shù)十倍。若只計(jì)算個(gè)別組件，計(jì)算效率還能提升。

MCNTC給出79盒燃料組件中各棒的功率分布，選取1／3堆芯，MCNTC與MCNP的計(jì)算結(jié)果對(duì)比如圖4所示。

圖4 MCNTC與MCNP計(jì)算的組件功率分布結(jié)果對(duì)比Fig.4 Comparison of power distributions for MCNTC and MCNP calculation

圖4中：MAX為棒相對(duì)功率最大偏差，按式（11）計(jì)算；AVE為棒功率代數(shù)平均偏差，按式（12）計(jì)算；RMS為棒功率幾何平均偏差，按式（13）計(jì)算；PDIF為最高功率棒偏差，按式（14）計(jì)算。

對(duì)計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)，MCNTC與MCNP計(jì)算結(jié)果的代數(shù)平均偏差小于0.8%，棒相對(duì)功率最大偏差小于1.9%，對(duì)79盒組件的統(tǒng)計(jì)數(shù)據(jù)列于表1。

表1 MCNTC與MCNP計(jì)算結(jié)果比較Table 1 Comparison of calculation results for MCNTC and MCNP

圖5 組件內(nèi)的相對(duì)功率分布與相對(duì)偏差Fig.5 Relative power distribution and relative deviation of assembly

圖5示出采用MCNTC計(jì)算得到的2－1、6－29位置組件內(nèi)相對(duì)功率分布及與MCNP計(jì)算結(jié)果的棒功率相對(duì)偏差。由圖5可看出，組件功率峰通常出現(xiàn)在燃料的邊角位置，且靠近堆芯中央。同時(shí)，組件位置越遠(yuǎn)離堆芯，組件功率峰因子越高。

表2列出部分組件位置MCNTC與NASPIN［10］計(jì)算結(jié)果的對(duì)比。表2中給出的結(jié)果均為與參考解的比值，參考解由MCNP計(jì)算。由表2可知：對(duì)于AVE，NAS－PIN的計(jì)算結(jié)果較MCNTC的稍好；對(duì)于PDIF，MCNTC的計(jì)算結(jié)果普遍較NAS－PIN的好。

3.2 CEFR核反應(yīng)率分布驗(yàn)證

CEFR在物理啟動(dòng)實(shí)驗(yàn)過(guò)程中，進(jìn)行了核反應(yīng)率分布測(cè)量［11］。實(shí)驗(yàn)過(guò)程中發(fā)現(xiàn)，197Au的俘獲相對(duì)反應(yīng)率通過(guò)MCNP程序F4卡獲得的理論計(jì)算值與實(shí)驗(yàn)值在反射層區(qū)存在較大偏差。經(jīng)過(guò)理論分析認(rèn)為，造成該偏差的主要原因是197Au的輻射俘獲截面在5eV附近存在一很強(qiáng)的共振峰，處于該共振峰峰值處的中子自由程與探測(cè)片厚度相近，故必須考慮該空間自屏影響，而之前的理論計(jì)算中未考慮該影響，導(dǎo)致在反射層區(qū)域計(jì)算反應(yīng)率約是實(shí)驗(yàn)值的2倍。

表2 MCNTC與NAS－PIN計(jì)算結(jié)果對(duì)比Table 2 Comparison of calculation results for MCNTC and NAS－PIN

采用MCNTC進(jìn)行核反應(yīng)率分布計(jì)算。其中，燃料區(qū)輻照組件橫截面如圖6所示。在輻照組件輻照管中部放置一直徑為5mm、厚度為30μm的金箔。堆芯區(qū)輻照組件位置為2－2、3－3、4－3、5－3、6－4，反射層區(qū)域輻照組件位置為7－5、9－6、11－8。燃料區(qū)輻照組件抽樣108次，反射層區(qū)輻照組件抽樣2×108次。

不同測(cè)量位置下，探測(cè)片相對(duì)反應(yīng)率分布如圖7所示，實(shí)驗(yàn)值與計(jì)算值分別以2－2位置歸一。由圖7可見(jiàn)：計(jì)算值與實(shí)驗(yàn)值符合較好，最大相對(duì)偏差小于5%；197Au的俘獲反應(yīng)率在燃料區(qū)基本保持不變，而在反射層區(qū)明顯上升，這主要是因?yàn)樵搮^(qū)域有較多中子慢化到共振區(qū)能量，共振中子對(duì)反應(yīng)率貢獻(xiàn)很大。

圖6 燃料區(qū)輻照組件橫截面Fig.6 Layout of test assembly in fuel zone

圖7 不同位置的俘獲反應(yīng)率分布Fig.7 Distribution of capture reaction rate in different positions

3.3 CEFR金箔探測(cè)片共振自屏驗(yàn)證實(shí)驗(yàn)設(shè)計(jì)及結(jié)果

為進(jìn)一步驗(yàn)證金箔的共振自屏效應(yīng)，在2014年CEFR提升功率過(guò)程中，設(shè)計(jì)了探測(cè)片測(cè)量實(shí)驗(yàn)。該實(shí)驗(yàn)將7片直徑為9mm、不同厚度的金箔疊放在一起，置于輻照組件內(nèi)，輻照組件位于堆芯第2排反射層8－5位置。于2014年3月14日在CEFR上以0.1%額定功率輻照2h。7片金箔的厚度列于表3。

表3 金箔的厚度Table 3 Thickness of gold foil

輻照后將輻照組件拆解，利用高純鍺譜儀逐個(gè)測(cè)量探測(cè)片中198Au位于411.8keV處的特征峰活度，用于確定197Au的（n，γ）反應(yīng)率。實(shí)驗(yàn)測(cè)量的探測(cè)片相對(duì)比活度及MCNTC計(jì)算結(jié)果列于表4。表4中的結(jié)果均按照活度最小的4號(hào)探測(cè)片歸一。

表4 探測(cè)片活度實(shí)驗(yàn)值與計(jì)算值結(jié)果比較Table 4 Comparison of gold foil activities for calculation and test

由表4可知，MCNTC計(jì)算結(jié)果與實(shí)驗(yàn)值符合較好，相對(duì)偏差在10%以內(nèi)。同時(shí)，該實(shí)驗(yàn)驗(yàn)證了金箔在快堆反射層區(qū)域存在顯著的共振自屏效應(yīng)，探測(cè)片表層反應(yīng)率是內(nèi)部反應(yīng)率的2倍左右。

4 結(jié)論

1）與功率重構(gòu)方法相比，對(duì)于規(guī)則組件，本文方法所得結(jié)果與其相近，計(jì)算速度稍慢。但本文方法對(duì)組件幾何與材料的適應(yīng)能力更好，能方便計(jì)算各類復(fù)雜幾何組件。

2）與MCNP相比，在全堆精細(xì)功率計(jì)算上本方法可得到與之相近的結(jié)果，但計(jì)算速度提升明顯，且組件離堆芯越遠(yuǎn)，相對(duì)效率提升越明顯。

3）本文方法可用于某些特殊的應(yīng)用，如反射層中探測(cè)片反應(yīng)率計(jì)算。

［1］ WAGNER M R.Three－dimensional nodal and transport theory methods for hexagonal－Zgeometry［J］.Nuclear Science and Engineering，1989，103（4）：377－391.

［2］ 謝仲生，鄧力.中子輸運(yùn)理論數(shù)值計(jì)算方法［M］.西安：西北工業(yè)大學(xué)出版社，2005.

［3］ ASKEW J.A characteristics formulation of the neutron transport equation in complicated geometries，AEEW－M 1108［R］.UK：Atomic Energy Authorty，1972.

［4］ FORGET B.A three dimensional heterogeneous coarse mesh transport method for reactor calculations［D］.USA：School of Mechanical Engineering，Georgia Institute of Technology，2006.

［5］ 徐李，馬大園，施工，等.快堆三維六角形節(jié)塊法輸運(yùn)計(jì)算研究［J］.原子能科學(xué)技術(shù)，2013，47（2）：161－165.

XU Li，MA Dayuan，SHI Gong，et al.Research of 3－D hexagonal nodal transport method for fast reactor［J］.Atomic Energy Science and Technology，2013，47（2）：161－165（in Chinese）.

［6］ Team X－M C.MCNP：A general Monte Carlo N－particle transport code，Version 5，LA－UR－03－1987［R］.USA：Los Alamos National Laboratory，2003.

［7］ 許淑艷.蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用［M］.北京：原子能出版社，2006.

［8］ MACFARLAN R E.The NJOY nuclear data processing system［R］.USA：Los Alamos National Laboratory，1994.

［9］ 田和春.中國(guó)實(shí)驗(yàn)快堆堆芯設(shè)計(jì)說(shuō)明［R］.北京：中國(guó)原子能科學(xué)研究院，2002.

［10］曹攀.快堆組件內(nèi)精細(xì)功率分布的計(jì)算研究［D］.北京：中國(guó)原子能科學(xué)研究院，2012.

［11］范振東，陳曉亮，胡定勝，等.中國(guó)實(shí)驗(yàn)快堆核反應(yīng)率分布測(cè)量研究［J］.原子能科學(xué)技術(shù)，2013，47（增刊）：107－110.

FAN Zhendong，CHEN Xiaoliang，HU Dingsheng，et al.Experimental research of nuclear reaction rate in China Experimental Fast Reactor［J］.Atomic Energy Science and Technology，2013，47（Suppl.）：107－110（in Chinese）.

［12］曹攀，喻宏，徐李，等.快堆燃料組件內(nèi)精細(xì)功率分布的計(jì)算［J］.原子能科學(xué)技術(shù)，2013，47（增刊）：287－290.

CAO Pan，YU Hong，XU Li，et al.Research on pin power distribution in fuel subassembly of fast reactor［J］.Atomic Energy Science and Technology，2013，47（Suppl.）：287－290（in Chinese）.

Study on Nodal and Monte Carlo Coupling Method

WANG Xin－zhe，YU Hong，XU Li，HU Yun，YANG Xiao－yan
（Department of Fast Reactor Research and Design，China Institute of Atomic Energy，Beijing102413，China）

The conventional nodal method could not calculate the distribution of fine neutron flux density in assemblies.A coupling method，called incident angular flux nodal and Monte Carlo coupling（IFNMC）method，was put forward.The computer code was programmed to verify this method.The calculation results are consistent quite well with the reference data.The IFNMC method is both effective and accurate for fine neutron flux density and other parameter calculations in special assembly.

nodal method；Monte Carlo method；coupling；incident angular flux；fine neutron flux density

TL329.2

A

1000－6931（2015）02－0200－07

10.7538／yzk.2015.49.02.0200

2014－09－16；

2014－12－02

863計(jì)劃資助項(xiàng)目（2011AA050301）

王新哲（1989—），男，山東東營(yíng)人，研究實(shí)習(xí)員，反應(yīng)堆物理專業(yè)