核反應堆堆芯功率分布數值模擬

廉麗莉，劉凱旋

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

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

1 堆芯差分方程

1.1 邊界條件

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

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

圖1 反應堆堆芯燃料組件分布及邊界

表1 IAEA問題的反照率數值解[8]

在編制計算程序的過程中，對反照率邊界條件進行進一步簡化，綜合考慮補充邊界節點代數方程的方法和附加源項法[9]，假設在反射層中與堆芯緊鄰處添加節點，且該節塊內的系數與堆芯相同，同時假設其中子通量密度Φr與該緊鄰堆芯節塊內中子通量密度Φr關系為：

Φr=β*Φc

(1)

1.2 堆芯差分方程形式

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

(2)

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

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

圖2 1/4反應堆堆芯組件圖

式(3)中各項系數如下所示：

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號節塊，即整體的(1,1)節塊，考慮全反射邊界條件差分方程形式為(5)形式

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

(5)

各項系數只有ai, j發生變化為

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

(6)

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

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)

各項系數只有ai, j發生變化為

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

(8)

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

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

(9)

各項系數只有ai, j發生變化為

(10)

其余元胞各類節塊差分方程的推導方法可根據元胞1類推。

1.3 中子擴散方程方程組系數矩陣

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

圖3 系數矩陣構造

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

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

2 二群中子擴散方程數值算法

本文采用SOR內迭代算法，其基本算法流程如下：

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

(11)

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

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)||∞≤ε或者達到最大迭代次數，計算終止，返回Φ(k)。

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

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

模擬過程的源迭代收斂準則為：

(13)

3 計算結果及分析

3.1 計算所需數據確定

在中子擴散方程數值解法的程序中，需要利用的核參數主要有：

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

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

3)反照率β，選取方法和相應數值見1.1；

4)節塊長度Δl，作為劃分網格的主要依據；

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

6)三區平均裂變中子數ν，用于在源迭代中計算新的源項；

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

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

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

表2 快群核數據

表3 熱群核數據

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

對于反照率β值，其取值過程1.1中已經詳細介紹，在此不贅述。

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

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

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

3.2 計算結果及分析

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

表4 運行特征數據

快群、熱群內迭代以及源迭代的收斂精度分別取10-4，10-4，10-4；10-5，10-5，10-5；10-5，10-6，10-6。選取第6次數值模擬結果，如圖4與圖5所示。發現中子通量密度分布相對于實際峰因子很大，分析認為是因為本程序采取了各種簡化手段，如認為反射層中無中子慢化，將反射層物理參數與堆芯歸一化處理，在反射層應用添加結塊邊界處理方法，認為堆芯內無外中子源等。為簡化問題，計算核數據只是MOX燃料棒的核數據，沒有考慮水以及其他材料的慢化作用。

由圖6與圖7中的曲線可見快群中子內迭代次數隨源迭代次數的增加減小。這是因為隨著源迭代次數增加快中子分布越來越精確，所以每后一次內迭代的初始數據都要比前一次精確。

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

圖4 第6次快群中子分布數值模擬結果

圖5 第6次熱群中子分布數值模擬結果

圖6 第3次計算內迭代次數變化曲線

圖7 第4次計算內迭代次數變化曲線

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

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

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

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

4 結論

本文數值模擬了反應堆堆芯快群中子以及熱群中子分布情況，認為相對于實際情況峰因子過大是因為反照率邊界條件應用不當，不同富集度組件中子分布差異劇烈是因為截面數據選取不當。

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