核反應(yīng)堆堆芯功率分布數(shù)值模擬

廉麗莉，劉凱旋

(中國原子能科學(xué)研究院，北京 102413)

中子擴(kuò)散方程的求解最常見的算法是有限差分法[1]，即通過對中子輸運(yùn)方程的各種簡化處理得到擴(kuò)散方程，再通過內(nèi)迭代外迭代劃分網(wǎng)格等得到最終的收斂解。隨著時代的發(fā)展，近年來，研究者開發(fā)了許多其他的用于求解中子擴(kuò)散方程解法，如節(jié)塊法[2]，現(xiàn)在已成為壓水堆設(shè)計中最為常用的方法。同時，國內(nèi)外有學(xué)者提出對中子擴(kuò)散方程進(jìn)行改造擴(kuò)展，引入兩個分別與中子吸收和中子擴(kuò)散有關(guān)的修正因子以彌補(bǔ)經(jīng)典中子擴(kuò)散方程的不足，從而可以有效計算多變量下的中子截面以及不均勻結(jié)構(gòu)間的中子通量密度分布[3-4]。除此之外，國外又發(fā)展出了中子擴(kuò)散方程的線性擴(kuò)展方法，利用這種方程得到的結(jié)果與半解析法得到的基準(zhǔn)數(shù)據(jù)吻合的很好，從而證明了擴(kuò)散方程線性擴(kuò)展的可行性[5]。作為中子擴(kuò)散方程的經(jīng)典解法，有限差分法至今仍有著很重要的地位，本文采用有限差分法數(shù)值模擬堆芯中子通量分布。

1 堆芯差分方程

1.1 邊界條件

在堆芯中子通量密度計算中，需要應(yīng)用到兩類邊界條件。首先因?yàn)槭菍?/4堆芯進(jìn)行計算，這樣在對稱分界面部分需要應(yīng)用全反射邊界條件，此類邊界條件較簡單在此不贅述。另一類邊界條件即反照率邊界條件，而反照率本身帶有很多經(jīng)驗(yàn)因子，無法精確得到其解析解，Segev[6]和M.Itagaki[7]分別提出應(yīng)用有限差分和邊界元數(shù)值法直接數(shù)值地產(chǎn)生邊界條件。由于相關(guān)文獻(xiàn)的報道很少，本文直接采用應(yīng)用于節(jié)塊法的一種邊界條件，如表1所示。

表1中界面AB等按照圖1對應(yīng)，其中，BCA，BCB分別為利用第一種和第二種方法得到的反照率數(shù)值解。

圖1 反應(yīng)堆堆芯燃料組件分布及邊界

表1 IAEA問題的反照率數(shù)值解[8]

在編制計算程序的過程中，對反照率邊界條件進(jìn)行進(jìn)一步簡化，綜合考慮補(bǔ)充邊界節(jié)點(diǎn)代數(shù)方程的方法和附加源項(xiàng)法[9]，假設(shè)在反射層中與堆芯緊鄰處添加節(jié)點(diǎn)，且該節(jié)塊內(nèi)的系數(shù)與堆芯相同，同時假設(shè)其中子通量密度Φr與該緊鄰堆芯節(jié)塊內(nèi)中子通量密度Φr關(guān)系為：

Φr=β*Φc

(1)

1.2 堆芯差分方程形式

式(2)為堆芯中子擴(kuò)散方程：

(2)

對于1/4堆芯組件，其節(jié)塊中子擴(kuò)散方程需要分為13種(如圖2所示)，首先分析元胞1節(jié)塊擴(kuò)散方程形式及規(guī)律。需要說明的是，本文的計算以組件均勻化系數(shù)為準(zhǔn)。假設(shè)現(xiàn)在已經(jīng)獲得堆芯組件擴(kuò)散系數(shù)1×47矩陣X，以及堆芯組件的移出截面數(shù)據(jù)1×47矩陣Y(編號以圖2為準(zhǔn))。

ai-1, jΦi-1, j+ai, j-1Φi, j-1+ai, jΦi, j+

ai, j+1Φi, j+1+ai+1, jΦi+1, j=Si, j

(3)

式(3)為普遍情況下節(jié)塊的差分方程。對于元胞1，即1～8號組件第1層節(jié)塊，上表面需要考慮全反射邊界條件，第150個節(jié)塊即(1，150)號節(jié)塊，則還要考慮反照率邊界條件。

圖2 1/4反應(yīng)堆堆芯組件圖

式(3)中各項(xiàng)系數(shù)如下所示：

ai+1, j=-2Di+1, jDi, j/Δl2(Di+1, j+Di, j)

ai-1, j=-2Di-1, jDi, j/Δl2(Di-1, j+Di, j)

ai, j-1=-2Di, j-1Di, j/Δl2(Di, j-1+Di, j)

ai, j+1=-2Di, j+1Di, j/Δl2(Di, j+1+Di, j)

ai, j=ΣR,i, j-ai-1, j-ai, j-1-ai, j+1-ai+1, j

(4)

元胞1的1號節(jié)塊，即整體的(1,1)節(jié)塊，考慮全反射邊界條件差分方程形式為(5)形式

ai, jΦi, j+ai, j+1Φi, j+1+ai+1, jΦi+1, j=Si, j

(5)

各項(xiàng)系數(shù)只有ai, j發(fā)生變化為

ai, j=ΣR,i, j-ai, j+1-ai+1, j

(6)

對于元胞1的第2～149號節(jié)塊即總體(1,2)～(1,149)號節(jié)塊，考慮全反射邊界條件差分方程形式為

ai, j-1Φi, j-1+ai, jΦi, j+ai, j+1Φi, j+1+ai+1, jΦi+1, j=Si, j

(7)

各項(xiàng)系數(shù)只有ai, j發(fā)生變化為

ai, j=ΣR,i, j-ai, j-1-ai, j+1-ai+1, j

(8)

對于元胞1的第150號節(jié)塊即總體第(1,150)節(jié)塊，考慮全反射及反照率邊界條件差分方程形式為

ai, j-1Φi, j-1+ai, jΦi, j+ai+1, jΦi+1, j=Si, j

(9)

各項(xiàng)系數(shù)只有ai, j發(fā)生變化為

(10)

其余元胞各類節(jié)塊差分方程的推導(dǎo)方法可根據(jù)元胞1類推。

1.3 中子擴(kuò)散方程方程組系數(shù)矩陣

對此矩陣進(jìn)行稀疏處理，系數(shù)矩陣構(gòu)造如圖3所示，M代表非零數(shù)字，其余位置皆為0。圖3給出了系數(shù)矩陣的具體構(gòu)成，由圖可見每一個括號表示一類元胞，共有13類元胞。

圖3 系數(shù)矩陣構(gòu)造

圖3中，括號旁邊的標(biāo)注數(shù)字為此類元胞包含的行數(shù)。單向箭頭所在行即大部分節(jié)塊共用差分方程行，其標(biāo)注數(shù)字為可用此方程代表的結(jié)塊數(shù)目。無單向箭頭的行即只代表一個節(jié)塊。雙向箭頭位于每一類元胞的第一行，表示這類元胞方程有系數(shù)項(xiàng)之間空隙包含的0位數(shù)目。橢圓圈住數(shù)字就是說明在此類元胞第一行方程的第一個非零系數(shù)前有多少個0空位。

考慮到系數(shù)矩陣方程組數(shù)目較多且均為手工編程組建，不可避免的會有錯誤。在一切正確的前提下，該系數(shù)矩陣是一個嚴(yán)格對角占優(yōu)陣，在迭代算法中，嚴(yán)格對角占優(yōu)陣有很好的收斂特性。因此額外編寫檢查矩陣的對角占優(yōu)性程序以檢查系數(shù)矩陣對角占優(yōu)性。

2 二群中子擴(kuò)散方程數(shù)值算法

本文采用SOR內(nèi)迭代算法，其基本算法流程如下：

1)取初始點(diǎn)Φ(0)，松弛因子w，置k=0，精度要求ε和最大迭代次數(shù)N；其中最佳松弛因子

(11)

ρ(BJ)為Jacobi迭代陣譜半徑，利用改進(jìn)乘冪法計算。

2)利用迭代式計算Φ(k+1)；SOR迭代式為

Φ(k+1)=BwΦ(k)+fw

(12)

其中，Bw=(D-wL)-1[(1-w)D+wU]，fw=w(D-wL)-1b。

3)若||Φ(k)-Φ(k-1)||∞≤ε或者達(dá)到最大迭代次數(shù)，計算終止，返回Φ(k)。

4)否則，k=k+1，轉(zhuǎn)(2)。

很明顯，由于A是嚴(yán)格對角占優(yōu)陣，其Jacobi迭代一定收斂也就是說其Jacobi迭代陣譜半徑一定小于1，由式(11)可知w一定小于2，因此，SOR迭代陣譜半徑一定小于1，則必收斂。

模擬過程的源迭代收斂準(zhǔn)則為：

(13)

3 計算結(jié)果及分析

3.1 計算所需數(shù)據(jù)確定

在中子擴(kuò)散方程數(shù)值解法的程序中，需要利用的核參數(shù)主要有：

1)三區(qū)散射系數(shù)D，用于差分方程系數(shù)的計算；

2)三區(qū)移出截面ΣR，用于差分方程對角系數(shù)的計算；

3)反照率β，選取方法和相應(yīng)數(shù)值見1.1；

4)節(jié)塊長度Δl，作為劃分網(wǎng)格的主要依據(jù)；

5)三區(qū)裂變截面Σf，用于在源迭代中計算新的裂變源項(xiàng)；

6)三區(qū)平均裂變中子數(shù)ν，用于在源迭代中計算新的源項(xiàng)；

7)三區(qū)快群到熱群的轉(zhuǎn)移截面Σs,1→2，用于熱群中子擴(kuò)散方程差分方程組的解。

8)三區(qū)中子裂變能譜χ，用于二群擴(kuò)散方程差分方程組的解。

本文采用MOX七群核燃料核數(shù)據(jù)[10]，處理得到二群核數(shù)據(jù)見表2與表3。表格中數(shù)據(jù)從左往右依次為擴(kuò)散系數(shù)、移出截面、裂變截面。

表2 快群核數(shù)據(jù)

表3 熱群核數(shù)據(jù)

對于節(jié)塊長度即網(wǎng)格的邊長Δl，考慮到一般不能大于中子擴(kuò)散長度的0.5～1倍，中子擴(kuò)散長度約為1.8 cm，以及燃料組件邊長為21.504 cm，確定節(jié)塊邊長為1.075 2 cm。

對于反照率β值，其取值過程1.1中已經(jīng)詳細(xì)介紹，在此不贅述。

對于裂變能譜χ，具體數(shù)值為χ1=1-1.176 1e-7，χ2=1.176 1e-7。

對于轉(zhuǎn)移截面Σs,1→2，1～3區(qū)的Σs,1→2分別為2.614 20e-3，2.533 10e-3，2.474 90e-3，單位為cm-1。

對于平均裂變中子數(shù)ν，三種富集度下分別對其去平均值為2.863 6，2.877 3，2.883 5，進(jìn)一步簡化取ν的三區(qū)平均值為2.874 8。

3.2 計算結(jié)果及分析

表4給出了6次運(yùn)行時的運(yùn)算特征數(shù)據(jù)，分別為運(yùn)算時間t/s，源迭代次數(shù)i_k，快群熱群SOR最佳松弛因子w1、w2，快群熱群Jacobi迭代陣最大特征值lamat_1、lamat_2，快群熱群乘冪法迭代次數(shù)i_lamat_1、i_lamat_2。

表4 運(yùn)行特征數(shù)據(jù)

快群、熱群內(nèi)迭代以及源迭代的收斂精度分別取10-4，10-4，10-4；10-5，10-5，10-5；10-5，10-6，10-6。選取第6次數(shù)值模擬結(jié)果，如圖4與圖5所示。發(fā)現(xiàn)中子通量密度分布相對于實(shí)際峰因子很大，分析認(rèn)為是因?yàn)楸境绦虿扇×烁鞣N簡化手段，如認(rèn)為反射層中無中子慢化，將反射層物理參數(shù)與堆芯歸一化處理，在反射層應(yīng)用添加結(jié)塊邊界處理方法，認(rèn)為堆芯內(nèi)無外中子源等。為簡化問題，計算核數(shù)據(jù)只是MOX燃料棒的核數(shù)據(jù)，沒有考慮水以及其他材料的慢化作用。

由圖6與圖7中的曲線可見快群中子內(nèi)迭代次數(shù)隨源迭代次數(shù)的增加減小。這是因?yàn)殡S著源迭代次數(shù)增加快中子分布越來越精確，所以每后一次內(nèi)迭代的初始數(shù)據(jù)都要比前一次精確。

熱群中子內(nèi)迭代過程只需幾次源迭代就可以得到比較精確的分布，其內(nèi)迭代次數(shù)在首次源迭代時相比較快群中子小很多，且可以發(fā)現(xiàn)熱群中子分布基本上就是快群中子分布按照富集度增加或減少得到的，這一現(xiàn)象通過其surf圖表現(xiàn)的很明顯。

圖4 第6次快群中子分布數(shù)值模擬結(jié)果

圖5 第6次熱群中子分布數(shù)值模擬結(jié)果

圖6 第3次計算內(nèi)迭代次數(shù)變化曲線

圖7 第4次計算內(nèi)迭代次數(shù)變化曲線

熱群中子擴(kuò)散方程的迭代式右端項(xiàng)b2=(xita2/keff(i_k))*S_origin+xigema_s_12.*Fi_1。其中xigema_s_12.*Fi_1的貢獻(xiàn)要明顯大于(xita2/keff(i_k))*S_origin項(xiàng)。這是因?yàn)樵错?xiàng)多為快中子，體現(xiàn)在xita2值非常小，本程序取值為xita2=1.7161e-7，幾近為0。這又反過來說明了內(nèi)迭代的初始迭代次數(shù)一定是很小且收斂速度一定是很快的。如此也很好的說明了快熱中子surf圖中分布形狀幾近相同這一情況。正是考慮到了熱群中子收斂很快這一現(xiàn)象故將其收斂精度由0.000 01調(diào)整至0.000 001。

由圖8與圖9中的曲線可以看出，有效增殖因數(shù)基本收斂到0.43左右，可以獲得穩(wěn)定源分布。

圖8 第3次計算Keff收斂過程

圖9 第4次計算Keff收斂過程

4 結(jié)論

本文數(shù)值模擬了反應(yīng)堆堆芯快群中子以及熱群中子分布情況，認(rèn)為相對于實(shí)際情況峰因子過大是因?yàn)榉凑章蔬吔鐥l件應(yīng)用不當(dāng)，不同富集度組件中子分布差異劇烈是因?yàn)榻孛鏀?shù)據(jù)選取不當(dāng)。

通過編程獲得了有效增殖因數(shù)、Jacobi迭代陣譜半徑、快熱群內(nèi)迭代次數(shù)的收斂曲線，本數(shù)值計算程序具有較快的收斂速度，計算結(jié)果合理。