秦帥, 李云召, 賀清明, 曹良志, 吳宏春

(西安交通大學 核科學與技術(shù)學院,陜西 西安 710049)

兩步法[1]是核反應堆物理工程計算中常用的計算方法,其計算精度取決于輸運計算產(chǎn)生的均勻化群常數(shù)是否可靠。根據(jù)均勻化理論[2],均勻化前后應保證積分反應率守恒和泄漏率守恒,其中,考慮均勻化區(qū)域的各向異性散射是實現(xiàn)泄漏率守恒的關(guān)鍵條件。目前的兩步法計算中,一般使用2類方法:1)產(chǎn)生均勻化區(qū)域的高階散射截面,并用于堆芯輸運計算中;2)產(chǎn)生均勻化區(qū)域的擴散系數(shù),并用于堆芯擴散計算中。其中,擴散系數(shù)的精確計算一般需要用到高階散射截面[3]。因此,產(chǎn)生精確的高階散射截面是考慮均勻化區(qū)域各向異性散射的基礎(chǔ)。

蒙特卡羅輸運計算方法[4]幾何處理能力強,基于連續(xù)能量模擬,可以直接考慮共振效應,計算精度很高,已在均勻化群常數(shù)的產(chǎn)生中得到了一定應用[5-8]。但是,蒙特卡羅方法在計算高階中子通量矩時統(tǒng)計漲落很大,一般以中子標通量為權(quán)重計算平均散射角余弦,帶來計算誤差[6]。使用蒙特卡羅方法計算的高階散射矩陣,在快中子反應堆問題中引起的有效增殖因子計算偏差達到500×10-5～600×10-5[9]。

本文提出一種基于中子均方位移守恒的平均散射角余弦計算方法,稱為中子均方位移守恒法。因為計算平均散射角余弦的主要目的是保證均勻化區(qū)域泄漏率守恒,而中子的泄漏直接取決于中子在均勻化區(qū)域內(nèi)的均方位移。該方法在蒙特卡羅輸運計算過程中對中子均方位移進行統(tǒng)計,在得到中子均方位移后,計算得到使其守恒的平均散射角余弦。基于西安交通大學核工程計算物理實驗室自主開發(fā)的蒙特卡羅粒子輸運計算程序NECP-MCX[10]對該方法進行了實現(xiàn),并使用各向異性較強的快中子反應堆問題對該方法進行了數(shù)值測試,驗證了其計算精度。

1 中子均方位移守恒法理論模型

1.1 高階散射截面與平均散射角余弦

高階散射截面是將散射角的分布使用勒讓德函數(shù)進行展開的系數(shù)矩,其定義為:

(1)

式中:r為空間變量;μ為散射角余弦;g′為入射能群;g為出射能群;Σs, g′→g(r,μ)為空間r處,以散射角余弦為μ,由g′群散射至g群的散射截面;Σs,g′→g(r)為空間r處,由g′群散射至g群的散射截面;Σs,n,g′→g(r)為空間r處,由g′群散射至g群的n階散射截面;fn,g′→g(r)為對應能群下散射角余弦分布函數(shù)的n階矩;Pn(μ)為n階勒讓德函數(shù);N為展開階數(shù)。

散射角余弦分布函數(shù)的n階矩為:

(2)

式中:E′和E分別為中子入射能量和出射能量;Eg′和Eg為能群g′和g的能量下界;φn(r,E′)為空間r處能量為E′的n階中子通量矩;fn(r,E′,E)為空間r處,能量由E′散射至E的散射角余弦分布函數(shù)n階矩,其計算公式為:

(3)

式中:f(r,E′,E,μ)為空間r處,中子以散射角余弦μ,能量由E′轉(zhuǎn)移至E的概率。

根據(jù)式(1)～(3),可以得到各階散射截面。對于大多數(shù)反應堆問題,一般取展開階數(shù)為1,一階散射截面可以寫為:

(4)

(5)

式中φ1(r,E′)為空間r處能量為E′的一階中子通量矩。

由式 (5)可以看出,精確的平均散射角余弦計算應以一階中子通量矩為權(quán)重,計算其分布函數(shù)的一階矩,但在蒙特卡羅輸運計算過程中,不同飛行方向的中子會發(fā)生相互抵消,導致一階中子通量矩的統(tǒng)計漲落很大。因此,在傳統(tǒng)計算中一般采取通量正比假設[11]:

φn(r,E′)=Cnφ(r,E′)

(6)

式中:φ(r,E′)為空間r處能量為E′的中子標通量;Cn為n階對應的常數(shù)。

根據(jù)式 (6),可以直接用中子標通量代替式 (5)中的一階中子通量矩,避免統(tǒng)計漲落過大的問題。但中子通量和一階中子通量矩并不滿足通量正比假設,從而引入了計算誤差。

1.2 中子均方位移與平均散射角余弦

計算高階散射截面和平均散射角余弦的主要目的是考慮中子的各向異性散射過程,最終影響中子在均勻化區(qū)域內(nèi)的直線飛行距離。因此,可以根據(jù)中子均方位移守恒,計算中子的平均散射角余弦,從而避免一階中子通量矩的計算。

1.2.1 中子均方位移計算[12]

首先,考慮單能中子在無限大均勻介質(zhì)內(nèi)飛行的情況。將中子由當前位置移動到下個碰撞點的過程稱為中子的一次飛行,那么中子每次飛行長度l的概率分布函數(shù)p(l)滿足:

p(l)dl=Σte-Σtldl

(7)

式中Σt為均勻介質(zhì)的總截面。

中子在若干次飛行后的均方位移可以表示為:

(8)

根據(jù)向量運算,式(8)為:

(9)

根據(jù)式 (7),式(9)中右端第1項可以寫為:

(10)

式中λ為平均自由程,λ=1/Σt。

式 (9)中右端第2項可以寫為[12]:

(11)

最終,中子在n次飛行后的均方位移可以表示為:

(12)

當飛行次數(shù)確定后,中子均方位移是平均散射角余弦的函數(shù)。一般情況下,該函數(shù)在[-1, 1]單調(diào)增長,即中子平均散射角余弦越大,中子在發(fā)生碰撞后越傾向于向前飛行,其均方位移就越大。

1.2.2 平均散射角余弦的計算

根據(jù)1.2.1節(jié)的推導,基于中子均方位移守恒,平均散射角余弦的計算方法為:

1)中子經(jīng)由源抽樣或碰撞后進入能群g,記錄其當前位置為rg,s,設置ng=0,wg=0;

2) 模擬中子輸運至rg,c處發(fā)生碰撞,更新該中子在當前能群的碰撞次數(shù)ng=ng+1;

3) 為處理隱俘獲,更新該中子的碰撞權(quán)重wg+1=wg+wi(wi為碰撞前的中子權(quán)重);

4) 檢查中子的出射能群,若中子能群不發(fā)生改變,則回到步驟2),若中子能群發(fā)生改變或被殺死,則更新計數(shù)Ng(ng)=Ng(ng)+Wg/ng及均方位移計數(shù)Dg(ng)=Dg(ng)+|rg,c-rg,s|2·Wg/ng;

5) 若中子仍存活或已模擬中子數(shù)未達到設定值,回到步驟1),否則退出。

全部模擬結(jié)束后,不同碰撞次數(shù)下的中子均方位移為:

(13)

式中Nc為用戶設定的最大飛行次數(shù)。

得到中子均方位移后,即可根據(jù)式 (12)搜索得到均勻化區(qū)域的平均散射角余弦。

需要注意,式 (12)基于無限大介質(zhì)得出的,但是,在兩步法計算中,經(jīng)常需要計算有限幾何區(qū)域的均勻化群常數(shù)。針對這種情況,本文采取的修正方法為基于該有限幾何區(qū)域的材料布置建立無限大問題,然后使用通量正比假設和本文方法分別計算無限大問題的平均散射角余弦,并計算修正因子:

(14)

(15)

另外需要注意,本文計算的平均散射角余弦僅對各群自散射截面進行修正,這是因為中子發(fā)生各向異性散射時,其能量損失相對較少,更易發(fā)生在自散射中,因此僅對自散射截面進行修正即可有效改善計算精度。最終的一階自散射截面計算公式為(忽略空間項):

(16)

2 中子均方位移守恒法數(shù)值測試

為了驗證本文方法,基于蒙特卡羅粒子輸運計算程序NECP-MCX開發(fā)了相應的中子均方位移統(tǒng)計功能。NECP-MCX為西安交通大學核工程計算物理實驗室自主開發(fā)的蒙特卡羅粒子輸運計算程序,已具備均勻化計算功能并在“華龍一號”啟動物理試驗中進行了驗證[8]。

2.1 一維問題

首先基于文獻[9]中的一維問題對本文方法進行測試,該問題來自文獻[13]中設計的低濃鈾鈉冷快堆,為簡化問題,對燃料區(qū)和反射層進行了打混。問題幾何如圖1所示,燃料區(qū)和反射層的核素組成如文獻[9]所示。

圖1 一維問題布置圖Fig.1 Layout of one-dimensional problem

本文使用NECP-MCX連續(xù)能量模式計算參考解,并在求解過程中分別統(tǒng)計燃料區(qū)和反射層區(qū)的33群均勻化群常數(shù)(能群結(jié)構(gòu)見表1),其中一階散射截面基于通量正比假設得到。計算時共投入1 050代,其中前50代為非活躍代,每代投入50萬個粒子。之后,使用中子均方位移守恒法分別計算了燃料和反射層的平均散射角余弦及修正因子,并對一階自散射截面進行修正。最后,分別使用修正和未修正的一階散射截面對該問題進行NECP-MCX多群計算。

實驗室利用Axis PTZ相機，通過圖像檢測Rovio的像素重心，然后由坐標變換函數(shù)轉(zhuǎn)化為攝像頭姿態(tài)(pan,tilt)。通過支持向量機算法間接實現(xiàn)(pan,tilt)和固定坐標系(x,y)的轉(zhuǎn)換。具體定位步驟如下：

表1 33群能群結(jié)構(gòu)Table 1 Structure for 33-group

表2中給出了一維問題有效增殖因子keff的計算結(jié)果,參考解為1.147 18±0.000 03?？梢钥闯?使用通量正比假設時,由于低估了系統(tǒng)的各向異性,多群計算的keff較大,使用本文方法對一階散射截面修正后可以有效改善多群計算結(jié)果。

表2 一維問題有效增殖因子計算結(jié)果

圖2為一維問題的中子通量分布比較。由于低估了各向異性,通量正比假設計算的反射層區(qū)域通量較大,而使用修正截面的計算結(jié)果則與參考解符合較好。圖3中給出了一維問題的中子通量相對偏差分布,修正方法將通量的最大相對偏差由9%降低為4%。

圖3 一維問題中子通量相對偏差分布Fig.3 Distribution of relative biases of neutron flux in one-dimensional problem

2.2 二維均勻堆芯問題

基于文獻[9]中的二維均勻堆芯問題對本文方法進行測試,為簡化問題,對燃料和反射層組件進行了打混,其核素組成與2.1小節(jié)相同。圖4中給出了本問題的堆芯布置,堆芯由5圈燃料組件和2圈反射層組件組成,其中組件對邊距為16.3 cm。

圖4 二維均勻堆芯問題布置圖Fig.4 Layout of two-dimensional homogeneous core problem

本文使用NECP-MCX連續(xù)能量模式計算參考解,并在求解過程中分別統(tǒng)計各位置處燃料組件和反射層組件的33群均勻化群常數(shù),其中一階散射截面基于通量正比假設得到。計算時共投入1 300代,其中前300代為非活躍代,每代投入50萬個粒子。然后,使用中子均方位移守恒法分別計算了燃料和反射層的平均散射角余弦及修正因子,并對一階自散射截面進行修正。最后,分別使用修正和未修正的一階散射截面對該問題進行NECP-MCX多群計算。

表3中給出了本問題的keff計算結(jié)果,參考解為1.089 78±0.000 07?？梢钥闯?使用修正截面可以將keff的偏差由0.007 96降低為-0.000 31。

表3 二維均勻堆芯問題有效增殖因子計算結(jié)果

圖5中給出了NECP-MCX計算得到的相對裂變反應率分布。圖6和圖7中分別給出了使用通量正比假設和修正截面計算得到的堆芯裂變反應率相對偏差分布?？梢钥闯?使用修正截面后,最大相對偏差由3.754%下降到-0.990%,相對均方根偏差由1.864%下降到0.612%。

圖5 二維均勻堆芯相對裂變反應率分布Fig.5 Distribution of relative fission reaction rates in two-dimensional homogeneous core

圖6 使用通量正比假設計算的二維均勻堆芯裂變反應率相對偏差分布Fig.6 Distribution of relative biases of fission reaction rates in two-dimensional homogeneous core calculated by using proportional flux assumption

圖7 使用修正截面計算的二維均勻堆芯裂變反應率相對偏差分布Fig.7 Distribution of relative biases of fission reaction rates in two-dimensional homogeneous core calculated by using corrected cross-sections

2.3 二維非均勻堆芯問題

基于經(jīng)合組織核能署(OECD/NEA)發(fā)布的鈉冷快堆基準題MET-1000問題[14],本文設計了二維非均勻堆芯問題。MET-1000為中等尺寸堆芯,堆芯內(nèi)裝載金屬燃料,組件對邊距為16.247 cm。為簡化問題,各類型組件使用均勻模型,組件內(nèi)不同材料核素組成及體積份額見文獻[14]。堆芯布置見圖8,堆芯內(nèi)插入控制棒,增加了系統(tǒng)的各向異性。

圖8 二維非均勻堆芯問題布置圖Fig.8 Layout of two-dimensional heterogeneous core problem

類似地,首先使用NECP-MCX連續(xù)能量模式計算參考解,并在求解過程中分別統(tǒng)計各位置處組件的24群均勻化群常數(shù)(能群結(jié)構(gòu)見表4),其中一階散射截面基于通量正比假設得到。計算時共投入1 300代,其中前300代為非活躍代,每代投入100萬個粒子。然后,使用中子均方位移守恒法分別計算了5種組件的平均散射角余弦及修正因子,并對一階自散射截面進行修正。最后,分別使用修正和未修正的一階散射截面對該問題進行NECP-MCX多群計算。

表4 24群能群結(jié)構(gòu)Table 4 Structure of 24-group

表5中給出了本問題的keff計算結(jié)果,參考解為0.950 61±0.000 04?？梢钥闯?使用修正截面可以將keff的偏差由0.007 17降低為0.002 66。

表5 二維非均勻堆芯問題有效增殖因子計算結(jié)果

圖9中給出了NECP-MCX計算得到的相對裂變反應率分布。

圖9 二維非均勻堆芯相對裂變反應率分布Fig.9 Distribution of relative fission reaction rates in two-dimensional heterogeneous core

圖10和圖11中分別給出了使用通量正比假設和修正截面計算得到的堆芯裂變反應率相對偏差分布?？梢钥闯?使用修正截面后,最大相對偏差由4.675%下降到0.920%,相對均方根偏差由2.444%下降到0.569%。

圖10 使用通量正比假設計算的二維非均勻堆芯裂變反應率相對偏差分布Fig.10 Distribution of relative biases of fission reaction rates in two-dimensional heterogeneous core calculated by using proportional flux assumption

圖11 使用修正截面計算的二維非均勻堆芯裂變反應率相對偏差分布Fig.11 Distribution of relative biases of fission reaction rates in two-dimensional heterogeneous core calculated by using corrected cross-sections

2.4 擴散計算結(jié)果

本文對修正方法計算得到的高階散射矩陣進行數(shù)值測試,也對基于該方法產(chǎn)生的擴散系數(shù)進行數(shù)值測試。基于Inflow近似[3]和斐克定律,可得[8]:

(17)

式中:φg為中子通量;Σt,g為總截面;Dg為擴散系數(shù)。

當蒙特卡羅輸運計算結(jié)束后,式 (17)中除擴散系數(shù)外各項均已知,該式成為了關(guān)于擴散系數(shù)的線性方程組,求解該線性方程組即可得到擴散系數(shù)。分別使用通量正比假設方法和修正方法產(chǎn)生式 (17)中的一階散射矩陣,計算得到多群擴散系數(shù),并基于該擴散系數(shù)對第2.3小節(jié)中的非均勻二維堆芯問題進行擴散計算,對基于本文方法產(chǎn)生的擴散系數(shù)進行數(shù)值測試。測試時計算條件設置與第2.3小節(jié)相同。

圖12和13為基于通量正比假設方法和修正方法的堆芯裂變反應率擴散計算相對偏差分布。

圖12 基于通量正比假設方法的堆芯裂變反應率擴散計算相對偏差分布Fig.12 Distribution of relative biases of fission reaction rates in core diffusion calculation based on proportional flux assumption

圖13 基于修正方法的堆芯裂變反應率擴散計算相對偏差分布Fig.13 Distribution of relative biases of fission reaction rates in core diffusion calculation based on correction method

可以看出,使用修正方法后,最大相對偏差由5.092%下降到1.214%,均方根偏差由2.609%下降到0.792%?？梢钥闯?修正方法計算得到的擴散系數(shù)更為精確。

3 結(jié)論

本文提出的基于中子均方位移守恒的平均散射角余弦計算方法,算例對該方法進行了數(shù)值測試,驗證了本文的計算精度。與傳統(tǒng)方法相比,中子均方位移守恒法的計算精度更高,可以有效降低強各向異性問題的計算偏差。

本文僅對一階自散射截面進行了修正,下一步可對完整的高階散射矩陣修正展開研究。