解析函數展開節塊法的數值穩定性分析

孫 偉，李 慶，王 侃

(1.清華大學 工程物理系，北京 100084；

2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610041)

方形幾何節塊法中，橫向積分是其中的關鍵技術，但在六角形幾何中直接采用橫向積分技術時橫向泄漏會出現奇異項，奇異項的出現給節塊中子注量率的求解帶來困難[1]。針對此問題，韓國KAIST的Cho等[2]提出了解析函數展開節塊(AFEN)法。AFEN法是在六角形幾何內不分離變量而直接利用嚴格滿足中子擴散方程的解析基函數，把節塊內的各群中子注量率近似展開。然后根據節塊間耦合條件，將節塊面平均凈流表示為節塊體平均中子注量率和面平均中子注量率的函數關系，代入節塊中子平衡方程形成全堆關于節塊體平均中子注量率和面平均中子注量率的方程組，分別用加速的高斯塞德爾方法和源迭代方法求解中子注量率和堆芯有效增殖因數keff。

研究發現利用AFEN法數值求解一些特殊問題時，如keff接近于某一節塊的無限介質增殖因數kinf(節塊在堆芯中實際凈泄漏接近于零)時，面平均凈流對中子注量率的關系系數出現震蕩，使得關于中子注量率方程組的系數矩陣出現病態，進而方程組的求解不收斂，無法得出中子注量率的解。方形幾何中，基于橫向積分技術的非線性迭代解析節塊法的這類不穩定性問題，Joo等[3]做了相關研究，對病態節塊采用線性近似或節塊展開法(NEM)代替解析節塊法求解；而六角形幾何中，針對中子注量率直接展開的AFEN法的這類不穩定性問題，尚未見相關研究，本文提出兩種解決方法來解決上述問題，即截斷近似方法和泰勒展開方法。

1 AFEN法數值不穩定性分析

1.1 AFEN法公式推導

對于反應堆內一個均勻的六角形節塊，其三維多群穩態中子擴散方程[1]可寫成如下形式：

(1)

其中：r為空間位置變量；V為節塊體積；Φ(r)=[Φ1(r)，Φ2(r)，…，ΦG(r)]T為G群中子注量率構成的向量。矩陣Λ(keff)的元素Λgg′為：

νΣfg′)/Dg

(2)

式(1)的解析解主要取決于矩陣Λ(keff)的特征值λm及相應特征向量um的數值類型。對于一般多群問題，λm既可能為實數也可能為復數，本文僅介紹λm均為實數的情況。

利用Λ(keff)的λm及um可將每群中子注量率展開為：

(3)

其中：U為矩陣Λ(keff)所有特征向量構成的矩陣；第m列對應于特征向量um；μgm為矩陣U的元素；Aml、Bml、Am0為展開系數；el為節塊展開時所選坐標軸上的單位向量，l取4時其具體表達式參見文獻[1]。

SN、CS為雙曲或三角函數：

(4)

利用凈流和中子注量率的關系可得節塊內面平均凈流與面平均中子注量率、體平均中子注量率的關系：

(5)

將式(5)代入式(6)可得全堆芯所有節塊關于面平均中子注量率、體平均中子注量率的方程組，可簡寫為式(7)，對式(7)的求解采用源迭代法，具體求解過程參見文獻[4]。

,2,…,G

(6)

MX=N

(7)

式中：系數矩陣M與群截面、堆芯幾何參數、組件不連續因子有關；X為堆芯所有面平均中子注量率構成的向量；N與體平均中子注量率有關。

1.2 AFEN法數值不穩定性原因探究

式(5)中系數矩陣F的每個元素都是雙曲或三角函數的表達式，與組件半寬度H、λm有關。表1列出矩陣元素F(G，1)(具體表達式詳見式(8)和(9)，其他元素的表達式在形式上與F(G，1)的相同)隨λm的變化情況(H取7.35)。可看出，當λm>0時，隨λm的減小，F(G，1)變化平緩，但λm減小到某一值時，F(G，1)出現畸變，在λm接近于0時F(G，1)變為無窮大；當λm<0時，在接近于0時F(G，1)也為無窮大，通過觀察F(G，1)的求解式可發現，當λm接近于0時，F(G，1)求解式各項分子、分母均接近于0，當兩個接近于0的數運算時會引入舍入誤差[5]。矩陣F元素的畸變會引起式(7)中系數矩陣M的元素也相差幾個數量級甚至為無窮大，導致M變為剛性矩陣，式(7)內迭代求解不收斂。

表1 F(G，1)值隨特征值的變化情況

當λm>0時，有：

(8)

當λm<0時，有：

F(G,1)=(H|λm|(-1+cos(H|λm|0.5))-

H2|λm|1.5sin(H|λm|0.5))/6(-2+(2+

H2|λm|)cos(H|λm|0.5))+|λm|0.5×

(-1+cos(H|λm|0.5))/6(-H|λm|0.5×

cos(H|λm|0.5)+sin(H|λm|0.5))+

H|λm|sin(H|λm|0.5)/3(-H|λm|0.5×

cos(H|λm|0.5)+sin(H|λm|0.5))+

|λm|0.5(H|λm|0.5cos(H|λm|0.5)×

(-1+H|λm|0.5cos(0.5H|λm|0.5))-

sin(H|λm|0.5))/2(-2+H2|λm|+

(2+H2|λm|)cos(H|λm|0.5)-

H|λm|0.5sin(H|λm|0.5))

(9)

系數矩陣F的元素值出現畸變，是由于矩陣Λ(keff)出現接近于0的特征值，本文僅給出2群截面時的公式推導解釋特征值接近于0的原因，2群時矩陣Λ(keff)的具體形式及截面間的轉換關系如下：

Λ(keff)=

(10)

Σr1=Σt1-Σ1→1;Σr2=Σt2-Σ2→2;

χ1=1;χ2=0;Σ2→1=0

(11)

其中，Σr1、Σr2為1、2群移除截面。由式(10)、(11)可知矩陣Λ(keff)行列式為：

(12)

當Λ(keff)出現接近于0的特征值時，Λ(keff)行列式接近于0，得：

·

(13)

又根據文獻[3]可知，2群截面時節塊的kinf為：

·

(14)

由式(13)、(14)可知，當keff接近于某一節塊的kinf時，矩陣Λ(keff)會出現接近于0的特征值，引起系數矩陣M的病態，出現不收斂情況。而kinf接近于keff的含義是節塊在堆芯中凈流為0，這多出現在堆芯中部，且堆芯含有多種富集度燃料，因為一般堆芯都有泄漏，凈流很難為0。

1.3 AFEN法的數值穩定性方法介紹

出現不收斂的直接原因是矩陣Λ(keff)的特征值很小時，M變為剛性矩陣，導致迭代不收斂，為解決這一問題，提出2種方法。

1) 截斷近似方法。從表1可知，在畸變前當特征值很小時，元素值基本不變，此時可將出現畸變的值截斷，直接用畸變前一刻的值代替。具體方法是當0<λm<10-6時，λm近似取為10-6，當-10-8<λm<0時，λm近似取為-10-8(由于λm在(-0.01，0.01)范圍內時元素值變化不大，且堆芯內出現特征值接近于0的節塊極少，因此截斷點的選擇基本不影響堆芯計算結果的精度，本文截斷點是參考F(G，1)變化趨勢給出的)。截斷近似的本質是對病態節塊采取近似常解析基函數展開。

2) 泰勒展開方法。當-10-8<λm<10-6時，對矩陣F元素求解式中的雙曲函數或三角函數進行泰勒展開，泰勒展開的目的是對求解式繼續化簡，避免出現分子、分母接近于0的項。表1列出5階和7階泰勒展開時元素F(G，1)的計算結果，可看出7階展開求出的F元素值更接近于真實值。泰勒展開的核心思想是高階多項式展開，與Joo等提出的對病態節塊直接采取NEM方法相比，7階泰勒展開保留原AFEN法的求解格式，在保證精度的同時更容易實現。

2 數值驗證

HANDF-D[4]是基于解析函數展開節塊法編制的可用于三維多群六角形幾何計算的程序，其正確性已得到多個基準題的驗證。為檢驗上述兩種穩定性方法的有效性，將其添加到HANDF-D程序中。驗證對象為改造后的三維VVER440基準題[1]，改造方法是減小VVER440材料1的兩群吸收截面(表2)，使keff接近于材料2的kinf(1.082 55)。將TRIVAC[6]的計算結果作為基準解，TRIVAC程序是由加拿大開發的可計算三維多群六角形幾何堆芯的中子學程序，采用了多種方法，包括有限差分方法、有限元方法等。

表2 VVER440材料1截面變化

在使用原HANDF-D程序計算改造后VVER440基準題(軸向分12層，每25 cm為1層)時，迭代中出現keff為1.082 56，接近于材料2的kinf(1.082 55)，矩陣Λ(keff)出現極小特征值1.429 2×10-7，迭代不收斂。表3列出采用兩種穩定性方法時keff的結果，可看出兩種方法計算均收斂且keff與TRIVAC的相比相對偏差僅為-0.02%。圖1示出功率分布的比較結果，功率分布最大相對偏差為2.03%，出現在燃料與反射層交界處，這與文獻[4]中指出的HANDF-D程序本身固有2.5%左右的功率偏差相當。從keff和功率分布結果可看出，兩種穩定性方法計算結果精確可靠，均能很好地解決解析函數展開節塊法計算不穩定性問題。

表3 改造后VVER440基準題keff計算結果

圖1 改造后VVER440基準題功率分布比較

3 結論

通過計算分析發現利用AFEN法求解六角形堆芯中子注量率，當keff接近于某一節塊kinf時，會使得此節塊的矩陣Λ(keff)出現極小特征值。利用此極小特征值求解全堆芯關于中子注量率方程組的系數矩陣時，由于引入舍入誤差使得系數矩陣變成剛性矩陣，導致迭代不收斂。針對此類問題，提出2種解決方法：截斷近似方法和泰勒展開方法。對改造后的VVER440基準題驗證表明，兩種方法均可有效解決不收斂問題，計算結果和參考解符合很好這一事實也充分說明這兩種方法有較高的精度。

參考文獻：

[1] 夏榜樣. 三維多群六角形節塊法及其應用研究[R]. 西安：西安交通大學，2006.

[2] CHO N Z, NOH J M. Analytic function expansion nodal method for hexagonal geometry[J]. Nucl Sci Eng, 1995, 121(3): 245-253.

[3] JOO H G, JIANG Guobing, DOWNAR T J. Stabilization techniques for the nonlinear analytic nodal method[J]. Nucl Sci Eng, 1998, 130(1): 47-59.

[4] 孫偉，倪東洋，李慶，等. 三維多群六角形幾何中子擴散程序開發[J]. 原子能科學技術，2013，47(10)：1 707-1 712.

SUN Wei, NI Dongyang, LI Qing, et al. Development of 3D multi-groups neutron diffusion code for hexagonal geometry[J]. Atomic Energy Science and Technology, 2013, 47(10): 1 707-1 712(in Chinese).

[5] 周鐵，徐樹方，張平文，等. 計算方法[M]. 北京：清華大學出版社，2006.

[6] HéBERT A, SEKKI D. A user guide for TRIVAC version4[R]. Montreal: Insitut de Genie Nucleaire Department Degenie Mecanique-Ecole Polytechnique de Montreal, 2009.