廉麗莉，劉凱旋

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

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

1 堆芯差分方程

1.1 邊界條件

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

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

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

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

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

Φr=β*Φc

(1)

1.2 堆芯差分方程形式

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

(2)

對于1/4堆芯組件，其節(jié)塊中子擴散方程需要分為13種(如圖2所示)，首先分析元胞1節(jié)塊擴散方程形式及規(guī)律。需要說明的是，本文的計算以組件均勻化系數(shù)為準。假設現(xiàn)在已經獲得堆芯組件擴散系數(shù)1×47矩陣X，以及堆芯組件的移出截面數(shù)據(jù)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)為普遍情況下節(jié)塊的差分方程。對于元胞1，即1～8號組件第1層節(jié)塊，上表面需要考慮全反射邊界條件，第150個節(jié)塊即(1，150)號節(jié)塊，則還要考慮反照率邊界條件。

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

式(3)中各項系數(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)

各項系數(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)

各項系數(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)

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

(10)

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

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

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

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

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

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

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

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

1)取初始點Φ(0)，松弛因子w，置k=0，精度要求ε和最大迭代次數(shù)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)||∞≤ε或者達到最大迭代次數(shù)，計算終止，返回Φ(k)。

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

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

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

(13)

3 計算結果及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對于轉移截面Σ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，進一步簡化取ν的三區(qū)平均值為2.874 8。

3.2 計算結果及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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