基于改進(jìn)多區(qū)delta-tracking方法的蒙特卡羅中子輸運(yùn)跟蹤與臨界計(jì)算驗(yàn)證

陳珍平，謝金森，郭 倩，張震宇，謝 超，于 濤,*

(1.南華大學(xué) 核科學(xué)技術(shù)學(xué)院，湖南 衡陽(yáng) 421001； 2.湖南省數(shù)字化反應(yīng)堆工程技術(shù)研究中心，湖南 衡陽(yáng) 421001； 3.南華大學(xué) 環(huán)境與安全工程學(xué)院，湖南 衡陽(yáng) 421001)

蒙特卡羅方法因其計(jì)算精度高、幾何適應(yīng)性強(qiáng)、收斂速度與問(wèn)題維數(shù)無(wú)關(guān)、具有優(yōu)良并行計(jì)算特性等優(yōu)勢(shì)而被廣泛應(yīng)用于粒子輸運(yùn)問(wèn)題的求解。因此，蒙特卡羅中子輸運(yùn)方法被國(guó)際反應(yīng)堆物理界認(rèn)為是新一代高精度堆芯物理分析方法的重要候選方法之一，也是近年來(lái)國(guó)際上的研究熱點(diǎn)[1-2]。目前，相對(duì)于確定論方法而言，蒙特卡羅中子輸運(yùn)方法的最大缺點(diǎn)在于計(jì)算效率較低，這也是阻礙蒙特卡羅中子輸運(yùn)方法進(jìn)行反應(yīng)堆工程化應(yīng)用的主要瓶頸之一。據(jù)相關(guān)研究表明[3-4]，蒙特卡羅中子輸運(yùn)計(jì)算中約30%～80%的計(jì)算時(shí)間用于中子輸運(yùn)跟蹤。因此，如何提高中子輸運(yùn)跟蹤的效率是提高蒙特卡羅中子輸運(yùn)計(jì)算效率的重要手段。為彌補(bǔ)傳統(tǒng)中子輸運(yùn)跟蹤方法帶來(lái)的潛在效率降低的缺陷，本文采用改進(jìn)多區(qū)delta-tracking方法引入虛截面來(lái)對(duì)模型進(jìn)行多區(qū)虛擬均勻化處理，進(jìn)而在中子輸運(yùn)跟蹤時(shí)則不考慮不同材料邊界穿越問(wèn)題，理論上可提高中子輸運(yùn)跟蹤效率。

1 蒙特卡羅中子輸運(yùn)方法

蒙特卡羅中子輸運(yùn)方法的本質(zhì)[5]就是模擬每個(gè)中子在給定模型系統(tǒng)中的隨機(jī)游動(dòng)歷史，即模擬中子從產(chǎn)生到消失的運(yùn)動(dòng)過(guò)程。1個(gè)中子在介質(zhì)系統(tǒng)中的游動(dòng)歷史可用1組狀態(tài)參數(shù)來(lái)描述，中子的狀態(tài)參數(shù)主要包括空間位置r、能量E和運(yùn)動(dòng)方向Ω，用狀態(tài)變量S=(r，E，Ω)表示。假設(shè)1個(gè)源中子初始狀態(tài)為S0=(r0，E0，Ω0)，中子在介質(zhì)材料中第m次碰撞后狀態(tài)為Sm=(rm，Em，Ωm)，其中rm為中子在第m次碰撞點(diǎn)的空間位置，Em為粒子在第m次碰撞后的能量，Ωm為粒子在第m次碰撞后的運(yùn)動(dòng)方向。因此，1個(gè)源中子在介質(zhì)系統(tǒng)中隨機(jī)游動(dòng)，經(jīng)過(guò)與介質(zhì)靶核的若干次碰撞后，其游動(dòng)歷史結(jié)束(如被吸收或逃出系統(tǒng))。因此，1個(gè)中子在系統(tǒng)中的游動(dòng)歷史可用1個(gè)狀態(tài)序列表示為：S0，S1，S2，…，SM-1，SM。因此，蒙特卡羅方法的本質(zhì)就是根據(jù)任意的中子狀態(tài)Sm確定其碰撞后的下一個(gè)狀態(tài)Sm+1，進(jìn)而確定每個(gè)中子的狀態(tài)序列(即游動(dòng)歷史)。如圖1所示，當(dāng)中子在系統(tǒng)中隨機(jī)游動(dòng)時(shí)，由狀態(tài)Sm=(rm，Em，Ωm)(m=1，2，…，M-1)經(jīng)過(guò)碰撞達(dá)到狀態(tài)Sm+1=(rm+1，Em+1，Ωm+1)時(shí)，需通過(guò)中子輸運(yùn)跟蹤確定其第m+1個(gè)狀態(tài)的空間位置rm+1、能量Em+1和運(yùn)動(dòng)方向Ωm+1。其中，中子輸運(yùn)跟蹤的1個(gè)重要任務(wù)就是確定中子從Sm狀態(tài)到達(dá)Sm+1狀態(tài)時(shí)的飛行距離。

圖1 傳統(tǒng)中子輸運(yùn)跟蹤方法Fig.1 Traditional neutron transport tracking method

假設(shè)中子以Sm狀態(tài)出發(fā)，以Em能量、沿Ωm方向飛行L長(zhǎng)度后，發(fā)生第m+1次碰撞時(shí)的位置為rm+1=rm+LΩm，其中L稱為中子兩次碰撞之間的輸運(yùn)長(zhǎng)度。由核反應(yīng)堆物理理論可知，輸運(yùn)長(zhǎng)度L服從的概率密度分布函數(shù)為：

f(L)=Σt(rm,Em)·

(1)

其中：Σt為介質(zhì)材料的宏觀總截面；l為中子的飛行徑跡長(zhǎng)度。對(duì)于式(1)，當(dāng)中子在均勻介質(zhì)中輸運(yùn)時(shí)，則輸運(yùn)長(zhǎng)度抽樣為：

(2)

其中，ξ是在0～1區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)。但是，當(dāng)中子在非均勻介質(zhì)系統(tǒng)中進(jìn)行輸運(yùn)時(shí)(圖1)，可能會(huì)依次穿過(guò)多層介質(zhì)才發(fā)生碰撞，根據(jù)指數(shù)分布的無(wú)記憶特點(diǎn)，此時(shí)需逐層考慮中子是否在該層介質(zhì)發(fā)生碰撞；如果不發(fā)生碰撞，則中子將沿當(dāng)前方向繼續(xù)飛行到下一層介質(zhì)邊界處，再判斷是否在該層介質(zhì)發(fā)生碰撞，直到中子在某一第I層介質(zhì)發(fā)生碰撞，此時(shí)輸運(yùn)長(zhǎng)度抽樣為：

(3)

其中：ΔLi為中子由rm出發(fā)、沿Ωm方向隨機(jī)游動(dòng)、在第i個(gè)介質(zhì)區(qū)域內(nèi)飛行的距離；Σt,I(Em)為中子發(fā)生碰撞反應(yīng)的第I層介質(zhì)的宏觀總截面。如圖1所示，在進(jìn)行非均勻介質(zhì)系統(tǒng)(如實(shí)際燃料組件或反應(yīng)堆堆芯)蒙特卡羅中子輸運(yùn)計(jì)算時(shí)，為確定式(3)的中子輸運(yùn)長(zhǎng)度，傳統(tǒng)中子輸運(yùn)跟蹤方法需根據(jù)中子射線方程與幾何模型的邊界面方程，聯(lián)立求解計(jì)算ΔLi。因此，在蒙特卡羅中子輸運(yùn)計(jì)算中，為確定式(3)的輸運(yùn)長(zhǎng)度，當(dāng)中子的平均自由程大于局部模型的宏觀尺寸時(shí)，頻繁大量的ΔLi距離計(jì)算是導(dǎo)致中子輸運(yùn)跟蹤效率降低的主要原因之一。

2 傳統(tǒng)單區(qū)delta-tracking中子輸運(yùn)跟蹤方法

為彌補(bǔ)在傳統(tǒng)中子輸運(yùn)跟蹤方法中，中子每次穿越材料邊界時(shí)都需要頻繁地進(jìn)行距離計(jì)算帶來(lái)效率降低的缺陷，Woodcock等[6]提出了delta-tracking方法來(lái)解決這一問(wèn)題，其基本原理是假設(shè)1個(gè)模型共包含N個(gè)幾何體，每個(gè)幾何體均被賦予了一種特定材料，整個(gè)模型中共有M種材料。傳統(tǒng)delta-tracking方法通過(guò)給每種材料的宏觀總截面添加1個(gè)虛截面，然后整個(gè)模型被虛擬均勻化為單一均勻介質(zhì)材料來(lái)進(jìn)行蒙特卡羅中子輸運(yùn)跟蹤(單區(qū)delta-tracking方法)，從而中子每步的輸運(yùn)長(zhǎng)度便可根據(jù)式(2)直接獲得，即使中子穿越材料邊界時(shí)，也不再需進(jìn)行邊界距離計(jì)算，從而節(jié)省了中子輸運(yùn)跟蹤時(shí)間，提高了中子輸運(yùn)計(jì)算效率。

圖2示出基于傳統(tǒng)單區(qū)delta-tracking的中子輸運(yùn)跟蹤方法。假設(shè)1個(gè)模型中包含M種介質(zhì)材料M1，M2，M3，…，MM，對(duì)于中子處于某一能量Ei，各材料在該能量Ei下的宏觀總截面為Σt，1，Σt，2，Σt，3，…，Σt，M。對(duì)所有材料對(duì)應(yīng)能量下的宏觀總截面進(jìn)行如下處理：

Σt,max=max{Σt,1,Σt,2,Σt,3,…,Σt,M}

(4)

Σt,max=Σt,1+Σ1,δ=

Σt,2+Σ2,δ=…=Σt,M+ΣM,δ

(5)

其中：Σt，max為當(dāng)前Ei能量下所有M種材料中宏觀總截面的最大值；Σ1，δ，Σ2，δ，Σ3，δ，…，ΣM，δ為虛截面，這些截面數(shù)據(jù)是人為定義和添加的。此時(shí)，整個(gè)模型被當(dāng)成一種均勻介質(zhì)來(lái)處理，可用式(2)直接進(jìn)行輸運(yùn)長(zhǎng)度抽樣，進(jìn)而避免了多次邊界距離計(jì)算，提高了中子輸運(yùn)跟蹤效率。

圖2 傳統(tǒng)單區(qū)delta-tracking中子輸運(yùn)跟蹤方法Fig.2 Neutron transport tracking based on traditional single-regional delta-tracking method

此時(shí)根據(jù)式(2)，中子每步的輸運(yùn)長(zhǎng)度抽樣為：

(6)

由于給每種材料的宏觀總截面引入了虛截面，因此式(6)中抽樣采用的Σt，max是包含虛截面部分的。因此為排除因虛截面的引入而帶來(lái)的與真實(shí)中子輸運(yùn)物理問(wèn)題的偏差，需將中子物理反應(yīng)分為真反應(yīng)和虛反應(yīng)(或假反應(yīng))，也就是按照式(6)進(jìn)行輸運(yùn)長(zhǎng)度抽樣并進(jìn)行中子輸運(yùn)跟蹤時(shí)，需采用拒絕抽樣來(lái)判斷真假中子物理反應(yīng)，其判斷條件如下。

1) 當(dāng)中子按式(6)輸運(yùn)長(zhǎng)度隨機(jī)游動(dòng)1步后，到達(dá)第i種介質(zhì)材料，然后判斷是否滿足條件：ξ1<Σt,i/Σt,max(Σt,i為中子當(dāng)前碰撞點(diǎn)處材料的宏觀總截面，ξ1為抽樣在0～1之間的隨機(jī)數(shù))；如果滿足該條件，則中子在該點(diǎn)發(fā)生真實(shí)物理反應(yīng)(真反應(yīng))，從而進(jìn)行后續(xù)的反應(yīng)核素、反應(yīng)類型、中子反應(yīng)后狀態(tài)進(jìn)行抽樣，然后進(jìn)行下一步輸運(yùn)。

2) 反之，如果不滿足上述判斷條件，則中子在該點(diǎn)發(fā)生虛反應(yīng)或假反應(yīng)，此時(shí)中子的能量和運(yùn)動(dòng)方向不發(fā)生變化，繼續(xù)以式(6)抽樣新的輸運(yùn)長(zhǎng)度(此時(shí)需抽樣新的隨機(jī)數(shù))，并以原來(lái)的能量沿原來(lái)運(yùn)動(dòng)方向向前飛行；然后返回步驟1重新進(jìn)行判斷，直到中子發(fā)生真實(shí)物理反應(yīng)、被吸收或離開(kāi)系統(tǒng)為止。

3 改進(jìn)多區(qū)delta-tracking中子輸運(yùn)跟蹤方法

3.1 基本原理

在傳統(tǒng)單區(qū)delta-tracking中子輸運(yùn)跟蹤方法中，對(duì)一個(gè)含有M種材料的模型，通過(guò)引入唯一一個(gè)當(dāng)前中子能量下的最大宏觀總截面Σt，max來(lái)對(duì)整個(gè)模型進(jìn)行單區(qū)虛擬均勻化處理，將非均勻介質(zhì)轉(zhuǎn)化為單一均勻介質(zhì)，此時(shí)將輸運(yùn)長(zhǎng)度抽樣公式由式(3)簡(jiǎn)化為式(6)以提高中子輸運(yùn)跟蹤效率。然而研究發(fā)現(xiàn)，當(dāng)模型局部區(qū)域存在強(qiáng)中子吸收體(或Σt，max較大的材料)時(shí)，此時(shí)模型部分材料區(qū)域(除強(qiáng)中子吸收體外的區(qū)域)滿足Σt,i≤Σt,max，則有Σt,i/Σt,max≈0，這時(shí)中子輸運(yùn)跟蹤過(guò)程進(jìn)行拒絕抽樣時(shí)會(huì)出現(xiàn)大量虛反應(yīng)，這一現(xiàn)象反而會(huì)導(dǎo)致delta-tracking跟蹤效率降低[7-8]。為彌補(bǔ)這一缺陷，本文提出了基于改進(jìn)多區(qū)delta-tracking的中子輸運(yùn)跟蹤方法。

通過(guò)對(duì)傳統(tǒng)單區(qū)delta-tracking方法分析，導(dǎo)致其效率降低的潛在因素在于，由于強(qiáng)中子吸收體(宏觀總截面較大的材料)的存在，導(dǎo)致強(qiáng)中子吸收體的Σt，max與其他材料的Σt，i相比量級(jí)大太多。為解決該問(wèn)題，采用改進(jìn)多區(qū)delta-tracking方法，根據(jù)模型中各區(qū)域材料對(duì)中子的吸收特性，對(duì)整個(gè)模型進(jìn)行區(qū)域分解，將模型劃分為多個(gè)虛擬均勻化區(qū)域(子區(qū)域)，在每個(gè)虛擬均勻化子區(qū)域內(nèi)再采用單區(qū)delta-tracking方法跟蹤。在進(jìn)行區(qū)域劃分時(shí)，直接利用構(gòu)成模型的幾何曲面進(jìn)行劃分，例如以燃料組件的邊界面(圖3)進(jìn)行虛擬劃分可得到多個(gè)子區(qū)域。此時(shí)，中子在不同子區(qū)域進(jìn)行delta-tracking跟蹤時(shí)，進(jìn)行輸運(yùn)長(zhǎng)度抽樣及拒絕抽樣時(shí)采用的Σt，max僅是當(dāng)前虛擬均勻化子區(qū)域內(nèi)包含材料的最大宏觀總截面，而每個(gè)子區(qū)域delta-tracking跟蹤時(shí)采用的Σt，max互不相同且互不影響。

a——組件模型；b——單區(qū)delta-tracking跟蹤；c——多區(qū)delta-tracking跟蹤圖3 改進(jìn)多區(qū)delta-tracking中子輸運(yùn)跟蹤方法Fig.3 Neutron transport tracking based on improved multi-regional delta-tracking method

以圖3a所示模型為例，其包含4個(gè)燃料組件，假設(shè)整個(gè)模型包含M種材料，其中1、2、3、4號(hào)組件各包含M1、M2、M3、M4種材料。在進(jìn)行delta-tracking中子輸運(yùn)跟蹤時(shí)，利用式(6)進(jìn)行輸運(yùn)長(zhǎng)度抽樣及拒絕抽樣時(shí)，對(duì)傳統(tǒng)單區(qū)delta-tracking跟蹤(圖3b)，其最大宏觀總截面為：

(7)

對(duì)于改進(jìn)多區(qū)delta-tracking跟蹤(圖3c)，假設(shè)將整個(gè)模型以組件為單位劃分為4個(gè)虛擬均勻化子區(qū)域進(jìn)行中子輸運(yùn)跟蹤，跟蹤時(shí)各子區(qū)域的最大宏觀總截面為：

(8)

通過(guò)分析式(7)、(8)可知，將單區(qū)改進(jìn)為多區(qū)跟蹤時(shí)，可減小Σt，max與Σt，i之間的量級(jí)差異。因此，通過(guò)將傳統(tǒng)delta-tracking方法的單區(qū)全局跟蹤方式改進(jìn)為基于多區(qū)的局部跟蹤方式，可有效避免因某個(gè)或幾個(gè)子區(qū)域存在強(qiáng)中子吸收體而導(dǎo)致其他區(qū)域中子輸運(yùn)跟蹤效率降低的現(xiàn)象。

3.2 無(wú)偏性證明

在進(jìn)行delta-tracking中子輸運(yùn)跟蹤時(shí)，當(dāng)中子穿越材料邊界時(shí)不再進(jìn)行式(3)的邊界距離計(jì)算，而此時(shí)的跟蹤過(guò)程與傳統(tǒng)蒙特卡羅中子輸運(yùn)跟蹤方法完全不同。但從概率論和統(tǒng)計(jì)學(xué)上來(lái)講，中子經(jīng)delta-tracking輸運(yùn)跟蹤得到的最終統(tǒng)計(jì)結(jié)果與采用傳統(tǒng)中子輸運(yùn)跟蹤方法所得到的結(jié)果是一致的，即delta-tracking方法對(duì)中子輸運(yùn)計(jì)算結(jié)果的統(tǒng)計(jì)是無(wú)偏的。以下證明給出了delta-tracking方法對(duì)輸運(yùn)結(jié)果統(tǒng)計(jì)的無(wú)偏性。

假設(shè)1個(gè)給定能量的中子在介質(zhì)中穿過(guò)路程L仍未發(fā)生反應(yīng)的概率為P(L)。由式(5)可知Σt,max=Σt+Σδ，則對(duì)該式兩邊從0到L進(jìn)行積分有：

(9)

在中子輸運(yùn)1步時(shí)，中子可能經(jīng)過(guò)0次，1次，2次，…，n次虛反應(yīng)到達(dá)L處，則有：

(10)

其中，P(L|n)可由如下得到：

P(L|0)=e-Σt,maxL

(11)

(12)

(13)

(14)

因此，中子穿過(guò)L長(zhǎng)的路程仍未發(fā)生核反應(yīng)的概率P(L)為：

(15)

由式(15)可知，在delta-tracking跟蹤方法中，中子輸運(yùn)長(zhǎng)度仍服從負(fù)指數(shù)分布，與傳統(tǒng)的蒙特卡羅中子輸運(yùn)跟蹤過(guò)程中，中子輸運(yùn)長(zhǎng)度服從的概率密度分布函數(shù)(式(1))在本質(zhì)上是一致的。因此，delta-tracking跟蹤方法對(duì)中子輸運(yùn)結(jié)果統(tǒng)計(jì)是無(wú)偏的。

3.3 中子輸運(yùn)跟蹤流程

基于改進(jìn)多區(qū)delta-tracking方法的中子輸運(yùn)跟蹤流程如圖4所示，給出了中子完成1步隨機(jī)游動(dòng)過(guò)程的輸運(yùn)跟蹤流程。

4 數(shù)值驗(yàn)證

本文將改進(jìn)多區(qū)delta-tracking中子輸運(yùn)跟蹤方法在NEAL團(tuán)隊(duì)開(kāi)發(fā)的多功能輻射輸運(yùn)模擬仿真平臺(tái)MOSRT系統(tǒng)[9]中進(jìn)行了程序?qū)崿F(xiàn)。本文開(kāi)展數(shù)值驗(yàn)證時(shí)分為兩個(gè)層面：為驗(yàn)證本文方法及程序的正確性與有效性，從國(guó)際臨界安全基準(zhǔn)評(píng)價(jià)實(shí)驗(yàn)手冊(cè)中選擇40道基準(zhǔn)題模型(多層介質(zhì)幾何模型)開(kāi)展臨界計(jì)算驗(yàn)證；為進(jìn)一步驗(yàn)證本文方法及程序?qū)τ趶?fù)雜問(wèn)題的適應(yīng)性，選取3個(gè)全堆芯模型開(kāi)展臨界計(jì)算驗(yàn)證，以證明本文方法及程序?qū)τ趶?fù)雜問(wèn)題同樣有效。其中，本文在Windows 7系統(tǒng)下采用單核Intel Core i5-5200 2.20 GHz處理器開(kāi)展數(shù)值驗(yàn)證，臨界計(jì)算條件為：共計(jì)算3 000個(gè)循環(huán)代，跳過(guò)500個(gè)非活躍代，每代5 000個(gè)中子，以計(jì)算得到的有效增殖因數(shù)(keff)作為正確性驗(yàn)證，以完成臨界計(jì)算的計(jì)算時(shí)間作為效率提升的有效性驗(yàn)證。

圖4 基于改進(jìn)多區(qū)delta-tracking的中子輸運(yùn)跟蹤流程Fig.4 Neutron transport tracking procedure based on improved multi-regional delta-tracking method

4.1 基準(zhǔn)題模型臨界驗(yàn)證

為減小計(jì)算結(jié)果的統(tǒng)計(jì)漲落，保證計(jì)算的準(zhǔn)確性，蒙特卡羅臨界計(jì)算keff時(shí)通常采用徑跡長(zhǎng)度估計(jì)、碰撞估計(jì)和吸收估計(jì)3種方法進(jìn)行組合計(jì)算[10]。然而，基于delta-tracking方法進(jìn)行中子輸運(yùn)跟蹤時(shí)，無(wú)法統(tǒng)計(jì)中子在每個(gè)材料區(qū)域單獨(dú)的徑跡長(zhǎng)度，導(dǎo)致計(jì)算keff時(shí)無(wú)法采用徑跡長(zhǎng)度估計(jì)進(jìn)行組合計(jì)算。因此，根據(jù)keff的物理定義，本文采用直接估計(jì)法[11](新生一代的裂變中子數(shù)目與其直屬上一代的中子數(shù)目之比)代替徑跡長(zhǎng)度估計(jì)法，以保證keff計(jì)算的準(zhǔn)確性。

國(guó)際臨界安全基準(zhǔn)評(píng)價(jià)實(shí)驗(yàn)手冊(cè)(ICSBEP)[12]是由OECD-NEA組織編撰，用于校驗(yàn)臨界計(jì)算方法和程序正確性的實(shí)驗(yàn)基準(zhǔn)題庫(kù)。為驗(yàn)證本文方法及程序?qū)崿F(xiàn)的正確性，從手冊(cè)中選取了40道基準(zhǔn)題，其主要為多層平板、球殼、圓柱及鈾水棒柵等臨界實(shí)驗(yàn)裝置，且覆蓋不同燃料富集度、燃料類型和中子能譜。其中，根據(jù)235U的富集度分成高富集度、中富集度和低富集度鈾；燃料類型包括金屬燃料和溶液燃料；中子能譜包括快譜和熱譜。因此，40道基準(zhǔn)題根據(jù)鈾富集度、燃料類型及中子能譜分為4類：1) hmf系列，高富集度、金屬燃料、快譜基準(zhǔn)題；2) hst系列，高富集度、溶液燃料、熱譜基準(zhǔn)題；3) imf系列，中富集度、金屬燃料、快譜基準(zhǔn)題；4) lst系列，低富集度、溶液燃料、熱譜基準(zhǔn)題。

表1列出基于多區(qū)delta-tracking(MDT)方法、單區(qū)delta-tracking(SDT)方法和傳統(tǒng)中子輸運(yùn)跟蹤方法(TTM)對(duì)40道基準(zhǔn)題進(jìn)行蒙特卡羅臨界計(jì)算得到的keff結(jié)果，并以MCNP程序計(jì)算結(jié)果作為基準(zhǔn)校驗(yàn)，同時(shí)給出實(shí)驗(yàn)值作為參考。從表1可知，基于MDT和SDT方法得到的keff計(jì)算結(jié)果與TTM結(jié)果的最大偏差分別為45 pcm和121 pcm，說(shuō)明MDT、SDT方法與TTM符合較好；同時(shí)，以MCNP結(jié)果作為基準(zhǔn)校驗(yàn)，MDT、SDT和TTM 3種方法得到的keff結(jié)果與MCNP結(jié)果的最大偏差分別為53、69和64 pcm，且3種方法的計(jì)算結(jié)果均落在MCNP結(jié)果的標(biāo)準(zhǔn)差內(nèi)，說(shuō)明在相同的臨界計(jì)算條件下3種方法與MCNP程序符合一致。因此，通過(guò)臨界基準(zhǔn)題的計(jì)算分析以及與MCNP程序的基準(zhǔn)校驗(yàn)，證明了本文方法及程序的正確性。

表2列出基于SDT和MDT方法進(jìn)行蒙特卡羅臨界計(jì)算時(shí)得到的加速比。從表2可知，基于delta-tracking方法可有效提升蒙特卡羅臨界計(jì)算的效率，可得到最大加速比為1.8。同時(shí)，通過(guò)臨界基準(zhǔn)題的計(jì)算驗(yàn)證發(fā)現(xiàn)，MDT方法得到的加速比略大于SDT方法，說(shuō)明MDT方法的中子輸運(yùn)跟蹤效率優(yōu)于SDT方法，證明了本文方法及程序在蒙特卡羅臨界計(jì)算效率提升方面的有效性。

4.2 全堆芯模型臨界驗(yàn)證

為進(jìn)一步驗(yàn)證本文方法及程序在復(fù)雜問(wèn)題計(jì)算中的適應(yīng)性，采用全堆芯模型對(duì)其進(jìn)行正確性和有效性驗(yàn)證，本文采用OPR、HM和BEAVRS 3個(gè)不同復(fù)雜度的全堆芯壓水堆(PWR)模型開(kāi)展了蒙特卡羅臨界計(jì)算驗(yàn)證。其中，OPR模型[13]由177個(gè)5種不同富集度燃料組件組成，每個(gè)燃料組件為16×16的棒柵布置；HM模型[14]由241個(gè)相同的燃料組件組成，每個(gè)燃料組件為17×17的棒柵布置；BEAVRS模型[15]由193個(gè)多種類型的燃料組件組成，每個(gè)燃料組件為17×17的棒柵布置。具體堆芯布置和結(jié)構(gòu)組成如圖5所示。

針對(duì)3個(gè)全堆芯模型進(jìn)行了臨界計(jì)算keff正確性驗(yàn)證。表3列出采用3種中子輸運(yùn)跟蹤方法得到的keff結(jié)果，并以MCNP程序計(jì)算結(jié)果作為基準(zhǔn)校驗(yàn)。從表3可知，基于MDT和SDT方法得到的keff結(jié)果與TTM結(jié)果的最大偏差分別為15 pcm和47 pcm，說(shuō)明MDT、SDT方法與TTM符合較好。同時(shí)以MCNP結(jié)果作為基準(zhǔn)校驗(yàn)，MDT、SDT和TTM 3種方法得到的keff結(jié)果與MCNP結(jié)果的最大偏差分別為18、46和33 pcm，且3種方法的計(jì)算結(jié)果均落在MCNP結(jié)果的標(biāo)準(zhǔn)差內(nèi)，說(shuō)明3種方法與MCNP程序符合一致。因此，數(shù)值結(jié)果證明了本文方法及程序在復(fù)雜全堆芯模型實(shí)際應(yīng)用中的正確性。

表4列出對(duì)3個(gè)全堆芯模型進(jìn)行臨界計(jì)算keff的有效性驗(yàn)證，以證明本文方法在處理復(fù)雜全堆芯模型時(shí)在中子輸運(yùn)跟蹤效率提升方面的有效性。從表4可知，MDT方法的中子輸運(yùn)跟蹤效率明顯優(yōu)于TTM和SDT方法，其對(duì)3個(gè)模型臨界計(jì)算的平均加速比約為1.4、1.5和1.6。因此，數(shù)值結(jié)果證明了本文方法及程序在復(fù)雜全堆芯模型實(shí)際應(yīng)用中的有效性。

5 結(jié)論

針對(duì)傳統(tǒng)蒙特卡羅中子輸運(yùn)跟蹤方法存在的效率降低問(wèn)題，本文提出了改進(jìn)多區(qū)delta-tracking中子輸運(yùn)跟蹤方法以提高中子輸運(yùn)跟蹤效率。本文對(duì)改進(jìn)多區(qū)delta-tracking方法的原理及跟蹤流程進(jìn)行了詳細(xì)闡述，并利用國(guó)際臨界安全實(shí)驗(yàn)基準(zhǔn)題和典型全堆芯壓水堆模型開(kāi)展了蒙特卡羅臨界計(jì)算驗(yàn)證。數(shù)值結(jié)果表明，改進(jìn)多區(qū)delta-tracking方法的臨界計(jì)算結(jié)果與傳統(tǒng)中子輸運(yùn)跟蹤方法的計(jì)算結(jié)果符合一致，且中子輸運(yùn)跟蹤效率得到明顯提升，證明了本文方法及程序的正確性和有效性。

表1 國(guó)際臨界安全基準(zhǔn)題模型臨界計(jì)算正確性驗(yàn)證Table 1 Verification of accuracy for ICSBEP benchmark criticality calculation

注：括號(hào)內(nèi)數(shù)據(jù)為標(biāo)準(zhǔn)差

表2 國(guó)際臨界安全基準(zhǔn)題模型臨界計(jì)算有效性驗(yàn)證Table 2 Verification of efficiency for ICSBEP benchmark criticality calculation

a——OPR模型；b——HM模型；c——BEAVRS模型圖5 全堆芯壓水堆模型Fig.5 Whole-core PWR model

全堆芯模型keffMDT方法SDT方法TTMMCNPOPRHMBEAVRS1.006 10(0.000 85)1.001 28(0.000 59)1.005 05(0.000 72)1.005 95(0.001 05)1.001 58(0.000 68)1.005 37(0.000 84)1.006 23(0.000 91)1.001 16(0.000 57)1.004 90(0.000 72)1.006 24(0.000 82)1.001 12(0.000 64)1.005 23(0.000 68)

注：括號(hào)內(nèi)數(shù)據(jù)為標(biāo)準(zhǔn)差

表4 全堆芯壓水堆模型臨界計(jì)算有效性驗(yàn)證Table 4 Verification of efficiency for whole-core PWR criticality calculation

注：1) 當(dāng)前中子輸運(yùn)跟蹤方法與傳統(tǒng)中子輸運(yùn)跟蹤方法的中子跟蹤效率之比

2) 傳統(tǒng)中子輸運(yùn)跟蹤方法與當(dāng)前中子輸運(yùn)跟蹤方法的總計(jì)算時(shí)間之比