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

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

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

圖2示出基于傳統(tǒng)單區(qū)delta-tracking的中子輸運(yùn)跟蹤方法。假設(shè)1個模型中包含M種介質(zhì)材料M1，M2，M3，…，MM，對于中子處于某一能量Ei，各材料在該能量Ei下的宏觀總截面為Σt，1，Σt，2，Σt，3，…，Σt，M。對所有材料對應(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ù)是人為定義和添加的。此時，整個模型被當(dāng)成一種均勻介質(zhì)來處理，可用式(2)直接進(jìn)行輸運(yùn)長度抽樣，進(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

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

(6)

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

1) 當(dāng)中子按式(6)輸運(yùn)長度隨機(jī)游動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)，此時中子的能量和運(yùn)動方向不發(fā)生變化，繼續(xù)以式(6)抽樣新的輸運(yùn)長度(此時需抽樣新的隨機(jī)數(shù))，并以原來的能量沿原來運(yùn)動方向向前飛行；然后返回步驟1重新進(jìn)行判斷，直到中子發(fā)生真實(shí)物理反應(yīng)、被吸收或離開系統(tǒng)為止。

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

3.1 基本原理

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

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

(7)

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

(8)

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

3.2 無偏性證明

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

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

(9)

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

(10)

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

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

(11)

(12)

(13)

(14)

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

(15)

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

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

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

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

本文將改進(jìn)多區(qū)delta-tracking中子輸運(yùn)跟蹤方法在NEAL團(tuán)隊(duì)開發(fā)的多功能輻射輸運(yùn)模擬仿真平臺MOSRT系統(tǒng)[9]中進(jìn)行了程序?qū)崿F(xiàn)。本文開展數(shù)值驗(yàn)證時分為兩個層面：為驗(yàn)證本文方法及程序的正確性與有效性，從國際臨界安全基準(zhǔn)評價實(shí)驗(yàn)手冊中選擇40道基準(zhǔn)題模型(多層介質(zhì)幾何模型)開展臨界計(jì)算驗(yàn)證；為進(jìn)一步驗(yàn)證本文方法及程序?qū)τ趶?fù)雜問題的適應(yīng)性，選取3個全堆芯模型開展臨界計(jì)算驗(yàn)證，以證明本文方法及程序?qū)τ趶?fù)雜問題同樣有效。其中，本文在Windows 7系統(tǒng)下采用單核Intel Core i5-5200 2.20 GHz處理器開展數(shù)值驗(yàn)證，臨界計(jì)算條件為：共計(jì)算3 000個循環(huán)代，跳過500個非活躍代，每代5 000個中子，以計(jì)算得到的有效增殖因數(shù)(keff)作為正確性驗(yàn)證，以完成臨界計(jì)算的計(jì)算時間作為效率提升的有效性驗(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時通常采用徑跡長度估計(jì)、碰撞估計(jì)和吸收估計(jì)3種方法進(jìn)行組合計(jì)算[10]。然而，基于delta-tracking方法進(jìn)行中子輸運(yùn)跟蹤時，無法統(tǒng)計(jì)中子在每個材料區(qū)域單獨(dú)的徑跡長度，導(dǎo)致計(jì)算keff時無法采用徑跡長度估計(jì)進(jìn)行組合計(jì)算。因此，根據(jù)keff的物理定義，本文采用直接估計(jì)法[11](新生一代的裂變中子數(shù)目與其直屬上一代的中子數(shù)目之比)代替徑跡長度估計(jì)法，以保證keff計(jì)算的準(zhǔn)確性。

國際臨界安全基準(zhǔn)評價實(shí)驗(yàn)手冊(ICSBEP)[12]是由OECD-NEA組織編撰，用于校驗(yàn)臨界計(jì)算方法和程序正確性的實(shí)驗(yàn)基準(zhǔn)題庫。為驗(yàn)證本文方法及程序?qū)崿F(xiàn)的正確性，從手冊中選取了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)對40道基準(zhǔn)題進(jìn)行蒙特卡羅臨界計(jì)算得到的keff結(jié)果，并以MCNP程序計(jì)算結(jié)果作為基準(zhǔn)校驗(yàn)，同時給出實(shí)驗(yàn)值作為參考。從表1可知，基于MDT和SDT方法得到的keff計(jì)算結(jié)果與TTM結(jié)果的最大偏差分別為45 pcm和121 pcm，說明MDT、SDT方法與TTM符合較好；同時，以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)，說明在相同的臨界計(jì)算條件下3種方法與MCNP程序符合一致。因此，通過臨界基準(zhǔn)題的計(jì)算分析以及與MCNP程序的基準(zhǔn)校驗(yàn)，證明了本文方法及程序的正確性。

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

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

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

針對3個全堆芯模型進(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，說明MDT、SDT方法與TTM符合較好。同時以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)，說明3種方法與MCNP程序符合一致。因此，數(shù)值結(jié)果證明了本文方法及程序在復(fù)雜全堆芯模型實(shí)際應(yīng)用中的正確性。

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

5 結(jié)論

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

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

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

表2 國際臨界安全基準(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)

注：括號內(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ì)算時間之比