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

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

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

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

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

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

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

高階散射截面是將散射角的分布使用勒讓德函數(shù)進(jìn)行展開(kāi)的系數(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)為對(duì)應(yīng)能群下散射角余弦分布函數(shù)的n階矩;Pn(μ)為n階勒讓德函數(shù);N為展開(kāi)階數(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階矩,其計(jì)算公式為:

(3)

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

根據(jù)式(1)～(3),可以得到各階散射截面。對(duì)于大多數(shù)反應(yīng)堆問(wèn)題,一般取展開(kāi)階數(shù)為1,一階散射截面可以寫為:

(4)

(5)

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

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

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

(6)

式中:φ(r,E′)為空間r處能量為E′的中子標(biāo)通量;Cn為n階對(duì)應(yīng)的常數(shù)。

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

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

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

1.2.1 中子均方位移計(jì)算[12]

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

p(l)dl=Σte-Σtldl

(7)

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

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

(8)

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

(9)

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

(10)

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

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

(11)

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

(12)

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

1.2.2 平均散射角余弦的計(jì)算

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

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

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

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

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

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

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

(13)

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

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

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

(14)

(15)

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

(16)

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

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

2.1 一維問(wèn)題

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

圖1 一維問(wèn)題布置圖Fig.1 Layout of one-dimensional problem

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

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

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

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

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

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

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

2.2 二維均勻堆芯問(wèn)題

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

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

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

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

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

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

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

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

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

2.3 二維非均勻堆芯問(wèn)題

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

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

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

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

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

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

圖9中給出了NECP-MCX計(jì)算得到的相對(duì)裂變反應(yīng)率分布。

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

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

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

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

2.4 擴(kuò)散計(jì)算結(jié)果

本文對(duì)修正方法計(jì)算得到的高階散射矩陣進(jìn)行數(shù)值測(cè)試,也對(duì)基于該方法產(chǎn)生的擴(kuò)散系數(shù)進(jìn)行數(shù)值測(cè)試?；贗nflow近似[3]和斐克定律,可得[8]:

(17)

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

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

圖12和13為基于通量正比假設(shè)方法和修正方法的堆芯裂變反應(yīng)率擴(kuò)散計(jì)算相對(duì)偏差分布。

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

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

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

3 結(jié)論

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

本文僅對(duì)一階自散射截面進(jìn)行了修正,下一步可對(duì)完整的高階散射矩陣修正展開(kāi)研究。