曹 璐，劉宙宇，賀清明，曹良志

(1.西安交通大學 核科學與技術學院，陜西 西安 710049；2.西北核技術研究所，陜西 西安 710024)

傳統(tǒng)共振計算方法包括等價理論[1]、子群方法[2]和超細群方法[3]，這些方法在處理復雜幾何燃料的共振計算時均存在各自的問題。等價理論存在幾何適應性問題，雖然通過三項有理近似可近似處理復雜幾何的逃脫概率，但計算精度不夠[4]；子群方法，理論上幾何適應性取決于子群固定源方程求解器的幾何處理能力，而隨著計算機性能的提升，特征線方法(MOC)廣泛應用于各確定論軟件中，因此子群方法可借用特征線方法的幾何適用性處理復雜幾何燃料，但由于子群方法計算中需要利用共振積分表，采用均勻的積分表計算精度不夠，而采用非均勻的共振積分表，需要提前針對幾種固定幾何的燃料進行共振積分表的計算，因此也存在幾何適應性問題；直接超細群方法僅用于小規(guī)模問題的計算，如用于作共振積分表的均勻問題、一維燃料棒和一維燃料板，這是由于在應用超細群方法時，需要計算區(qū)域之間的碰撞概率，對于一維燃料棒或一維燃料板問題，其區(qū)域之間的碰撞概率可采用解析解獲得，而對于復雜幾何燃料，碰撞概率的計算存在困難。國際上有研究[5]將特征線方法應用于超細群慢化方程的求解，以解決超細群方法的幾何適應性問題，但在應用該方法時，又需要提前作表用于加速，而作表又使得該方法只能適用于作表中采用的幾何針對高保真數(shù)值反應堆技術對計算精度和計算效率的要求，國際上對全堆芯高效、高精度共振計算方法開展了大量的研究[6-7]。NECP-X[8]是由西安交通大學核工程計算物理實驗室(NECP)自主開發(fā)的數(shù)值反應堆物理計算程序，已實現(xiàn)棒狀、板狀等幾何形式燃料的全堆芯高精度共振計算[9-10]，但對于非規(guī)則幾何燃料(如蜂窩煤狀燃料)，由于采用特征線方法[11]求解偽核素子群方程，因此計算效率很低。本研究提出一種基于等效幾何理論的共振計算方法用于處理各種幾何形式的燃料，并將該方法應用于NECP-X，以解決復雜幾何燃料的計算效率問題。

1 理論模型

1.1 基于等效幾何理論的共振計算方法

共振計算的目的是獲得共振能區(qū)的中子注量率，歸并獲得共振能量段的有效自屏截面，該過程如式(1)所示：

(1)

式中：下標x和g分別為反應道和能群編號；r、Ω和E分別為空間變量、角度變量和能量變量；ΔEg為能群的能量區(qū)間；σ和φ分別為微觀截面和中子注量率。

共振計算的難點在于能群內(nèi)中子通量密度的獲得，對于多區(qū)系統(tǒng)，中子通量密度可通過慢化方程求解：

Σt,f(E)φf(E)Vf=Pf→f(E)Sf(E)Vf+

(2)

式中：f為目標燃料區(qū)編號；f′為其他燃料區(qū)編號；m為慢化劑區(qū)編號；E為能量點，eV；F為所有燃料區(qū)集合；M為所有慢化劑區(qū)集合；Σt,f為燃料區(qū)的宏觀總截面，cm-1；φf為目標燃料區(qū)的中子通量密度，cm-2·s-1；V為區(qū)域體積，cm3；S為區(qū)域中子源強度，cm-3·s-1；Pf′→f為從其他燃料區(qū)產(chǎn)生的中子到目標燃料區(qū)第1次發(fā)生碰撞的概率；Pm→f為從慢化劑區(qū)產(chǎn)生的中子到目標燃料區(qū)第1次發(fā)生碰撞的概率。

基于互易公式，對慢化劑區(qū)源項采用窄共振近似，式(2)可寫為：

Σt,f(E)φf(E)=[1-Pf→M(E)]Sf(E)+

Pf→M(E)Σt,f(E)/E

(3)

式中：Pf→M(E)為目標燃料區(qū)產(chǎn)生的中子到慢化劑區(qū)第1次發(fā)生碰撞的概率。

由式(3)可知，對于不同幾何的燃料，如果能保證目標燃料區(qū)產(chǎn)生的中子在慢化劑中發(fā)生首次碰撞的概率Pf→M(E)相等，則能保證不同系統(tǒng)中的中子注量率相等，由式(1)可知，不同燃料的有效自屏截面將相同。

另外，根據(jù)丹可夫修正因子的定義可得：

(4)

式中：C為丹可夫修正因子；Pe,iso(E)為孤立系統(tǒng)下目標燃料區(qū)的逃脫概率。

所以Pf→m(E)可表示為：

Pf→m(E)=[1-C(E)]Pe,iso(E)

(5)

因此，如果保證不同幾何燃料中的(1-C(E))和Pe,iso(E)相等，則不同幾何燃料可獲得相同的有效自屏截面。

1.2 逃脫概率的計算

(6)

后續(xù)有研究采用貝爾因子aB對式(5)進行改進：

(7)

在輕水堆計算中，貝爾因子的選區(qū)范圍為1.1～1.4。

為進一步提高計算精度，Carlvik等[13]針對棒狀幾何燃料提出了兩項有理近似公式：

(8)

針對板狀幾何燃料，Roman等[14]同樣給出了兩項有理近似表達式：

(9)

但對于復雜幾何的燃料，難以獲得簡單直接的逃脫概率計算公式，而是通過求解慢化方程獲得。孤立系統(tǒng)中，根據(jù)中子守恒可得：

Σt,f(E)φf(E)Vf=Pf→f(E)Sf(E)Vf+

Pm→f(E)Sm(E)Vm

(10)

采用互易關系，式(10)可轉化為：

Σt,f(E)φf(E)=Pf→f(E)Sf(E)+

(11)

對于孤立系統(tǒng)，Pf→m(E)=Pe,iso(E)，所以式(10)可寫為：

(12)

通過上述推導可知，在孤立系統(tǒng)下可通過計算燃料區(qū)的中子注量率及式(12)獲得燃料區(qū)的逃脫概率。所以對于復雜幾何燃料可構造孤立系統(tǒng)問題，然后采用可處理復雜幾何的MOC求解該孤立問題用以獲得燃料區(qū)中子注量率進而獲得該幾何燃料的逃脫概率，再將其轉化為逃脫概率相同的一維圓柱(或平板)燃料模型。

以蜂窩煤狀燃料為例，基于逃脫概率守恒將蜂窩煤狀燃料等效為棒狀燃料的過程示于圖1。對于蜂窩煤狀燃料，采用MOC求解該燃料的孤立問題獲得該燃料的逃脫概率Pe,iso,1(E)，然后進行不同半徑下相同材料的棒狀燃料進行孤立問題的計算，獲得Pe,iso,2(E,r)，當Pe,iso,2(E,r)=Pe,iso,1(E)時所獲得的燃料半徑即為等效棒狀燃料的半徑。

圖1 基于等效幾何理論的蜂窩煤狀燃料的等效過程

1.3 燃料外圍結構材料的等效處理

燃料外圍的結構材料會影響燃料區(qū)的能譜，因此需要考慮結構材料的影響。對于復雜幾何燃料外圍的結構材料，通過保證燃料中產(chǎn)生的中子到結構材料的首次碰撞的概率守恒，在等效一維圓柱燃料外圍建立一層等效的結構材料。

在孤立系統(tǒng)中針對包殼區(qū)域建立中子平衡方程，可得：

Σt,clad(E)φclad(E)Vclad=Pf→clad(E)Sf(E)Vf+

(13)

式中：Σt,clad(E)為目標包殼區(qū)總截面；clad為包殼材料集。

基于互易關系(式(14))，式(13)可寫為式(15)。

Σt,clad(E)Pclad→f(E)Vclad=Pf→clad(E)Σt,f(E)Vf

Σt,clad(E)Pclad→m(E)Vclad=Pm→clad(E)Σt,m(E)Vm

(14)

Pf→clad(E)=(Σt,clad(E)φclad(E)Vclad-

(15)

從式(15)可看出，若要保證不同幾何下Pf→clad(E)守恒，則需要保證不同幾何下φclad(E)、Vclad/Vf、Pclad→m(E)Sm(E)/Σt,m、Pclad1→clad(E)·Sclad1(E)以及Sf(E)相等。對于原來的復雜幾何燃料的孤立問題，Pf→clad(E)可通過特征線方法獲得。

基于燃料區(qū)到包殼材料區(qū)碰撞概率守恒實現(xiàn)等效為棒狀燃料(蜂窩煤狀燃料)周圍包殼區(qū)的過程示于圖2。首先對圖2a目標燃料周圍的包殼采用MOC求解該燃料的孤立問題獲得該燃料到包殼的碰撞概率，然后在1.2節(jié)中已獲得的燃料棒半徑的基礎上，基于燃料到周圍結構材料碰撞概率守恒采用二分法進行等效結構材料半徑的搜索，獲得的等效結構材料幾何如圖2b所示。

圖2 基于等效幾何理論的蜂窩煤狀燃料的包殼區(qū)等效過程

1.4 丹可夫修正因子的處理

從1.1節(jié)可知，等效一維幾何燃料的丹可夫修正因子也應等于復雜幾何燃料在堆芯中的丹可夫因子。首先采用中子流方法獲得原復雜幾何燃料的丹可夫修正因子C(E)，然后搜索等效一維幾何燃料外圍的慢化劑尺寸，使C(pitch,E)=C(E)，C(pitch,E)為等效后幾何在慢化劑尺寸為pitch時的丹可夫修正因子。由上文可知，丹可夫修正因子是能量相關的，即灰體丹可夫修正因子，有研究[15]表明，采用黑體丹可夫修正因子C(Black)和灰體丹可夫修正因子C(E)時，1-C(Black)較1-C(E)偏差小，對最終有效自屏截面的計算影響很小，所以本文采用黑體丹可夫修正因子進行模擬計算。

基于燃料區(qū)丹可夫修正因子守恒實現(xiàn)等效為棒狀燃料(蜂窩煤狀燃料)外圍慢化劑的過程示于圖3。根據(jù)圖3a目標燃料全局計算獲得的丹克夫修正因子，在1.2、1.3節(jié)中獲得的燃料棒與周圍結構材料半徑的基礎上，基于丹可夫修正守恒采用二分法進行等效慢化劑半徑的搜索，獲得的等效結構材料幾何如圖3b所示。

圖3 基于等效幾何理論的蜂窩煤狀燃料的慢化劑等效過程

2 數(shù)值驗證

本文在NECP-X程序中實現(xiàn)基于等效幾何理論的復雜幾何燃料共振計算方法。為分析計算精度，對多種非棒狀柵元、組件問題進行計算分析。NECP-X采用基于ENDF/B-Ⅶ.0的拓展共振能群范圍的WIMS格式69群核數(shù)據(jù)庫，利用蒙特卡羅程序獲得參考解。NECP-X程序的計算條件為：1/4卦限內(nèi)有8個方位角和3個極角；特征線間距為0.01 cm；特征值收斂限為10-6；裂變率收斂限為10-5。為驗證本文方法的計算精度和計算效率，將本文方法計算結果與可處理復雜幾何燃料的傳統(tǒng)子群方法以及基于MOC的偽核素子群方法[10]的計算結果進行比較。

2.1 非棒狀燃料柵元問題

1)環(huán)板燃料柵元問題

環(huán)板燃料柵元問題燃料幾何取自于ATR反應堆，該環(huán)板燃料柵元問題的幾何以及平源區(qū)網(wǎng)格示于圖4，計算中省略包殼的建模，燃料采用3.1%二氧化鈾燃料，慢化劑采用輕水，環(huán)板燃料內(nèi)外徑分別為7.7 cm和7.750 8 cm，柵元長為6 cm，寬為1.6 cm，所有材料溫度均為600 K，采用全反射邊界條件，蒙特卡羅程序統(tǒng)計特征值結果為0.404 19±5.0×10-5。不同方法計算得到的235U和238U微觀吸收截面及其偏差示于圖5，不同共振計算方法的計算時間及特征值計算結果列于表1。

圖4 環(huán)板燃料柵元問題的幾何以及平源區(qū)網(wǎng)格

由圖5可見，基于MOC的偽核素子群方法計算的235U微觀吸收截面(σ)最大偏差發(fā)生在第27群，為1.69%，其他能群偏差均小于0.54%，238U微觀吸收截面最大偏差也發(fā)生在第27群，為1.56%，基于等效幾何理論的共振計算方法計算的235U微觀吸收截面偏差均小于0.9%，但相對于基于MOC的偽核素子群方法，238U微觀吸收截面偏差更大，最大偏差接近3.0%，傳統(tǒng)子群方法計算的235U微觀吸收截面最大偏差為5.67%，238U微觀吸收截面最大偏差為-1.75%，可看出，3種方法中基于MOC的偽核素子群方法精度最高，其次為基于等效幾何理論的共振計算方法，最后為傳統(tǒng)子群方法，但從表1的特征值可看出，3種方法的特征值計算結果相近，特征值偏差分別為84 pcm、62 pcm、52 pcm，上述截面及特征值結果表明，基于等效幾何理論的共振計算方法和基于MOC的偽核素子群方法在處理環(huán)板燃料共振計算中具有較好的計算精度，而從表1的共振計算時間可看出，采用基于等效幾何理論的共振計算方法的計算速度較基于MOC的偽核素子群方法快近80倍，較傳統(tǒng)子群方法快140.83倍，表明基于等效幾何理論的共振計算方法具有較高的計算效率。

圖5 環(huán)板燃料柵元問題中235U和238U吸收截面及偏差的計算結果

表1 環(huán)板燃料柵元問題中不同共振計算方法的計算時間及特征值計算結果

2)蜂窩煤狀燃料柵元問題

本文根據(jù)文獻[16]設計蜂窩煤狀燃料柵元問題，燃料采用3.1%二氧化鈾燃料，慢化劑采用輕水，所設計的蜂窩煤狀燃料柵元問題的幾何以及計算時采用的平源區(qū)網(wǎng)格示于圖6，蜂窩煤狀燃料內(nèi)部的每個水洞尺寸均為0.05 cm，共分布3圈19個水洞，3圈水洞圓心所在半徑分別為0.0、0.15、0.3 cm，蜂窩煤狀燃料外徑為0.4 cm，柵距為1.2 cm，所有材料溫度均為600 K，計算時采用全反射邊界條件。采用不同共振計算方法及蒙特卡羅程序計算235U和238U微觀吸收截面及其偏差，結果如圖7所示。由圖7可見，基于等效幾何理論的共振計算方法計算的235U的微觀吸收截面最大偏差為-0.48%，發(fā)生在第27群，238U的微觀吸收截面最大偏差為1.61%，發(fā)生在第24群；基于MOC的偽核素子群方法計算的235U的微觀吸收截面最大偏差為0.65%，發(fā)生在第27群，238U的微觀吸收截面最大偏差為1.88%，發(fā)生在第24群；傳統(tǒng)子群方法計算的235U的微觀吸收截面最大偏差為7.32%，發(fā)生在第25群，238U的微觀吸收截面最大偏差為1.31%，發(fā)生在第24群。特征值計算結果列于表2，蒙特卡羅程序統(tǒng)計結果為1.373 45±4.0×10-5，基于等效幾何理論的共振計算方法與基于MOC的偽核素子群方法計算的特征值偏差均小于50 pcm，而傳統(tǒng)子群方法計算的特征值偏差為-164 pcm；基于等效幾何理論的共振計算方法計算速度更快，相較于基于MOC的偽核素子群方法加速比為72.6，相較于傳統(tǒng)子群方法加速比為152.87。上述結果表明，基于等效幾何理論的共振計算方法在計算蜂窩煤狀燃料時有較好的計算精度和較高的計算效率。

圖6 蜂窩煤狀燃料柵元幾何及平源區(qū)網(wǎng)格

圖7 蜂窩煤狀燃料柵元問題中235U和238U吸收截面及偏差的計算結果

表2 蜂窩煤狀燃料柵元問題中不同共振計算方法的計算時間及特征值計算結果

2.2 非棒狀燃料組件問題

1)環(huán)板燃料組件問題

根據(jù)2.1節(jié)柵元問題建立環(huán)板燃料組件問題，共采用7個環(huán)板柵元，柵元幾何如圖8所示，所有材料的溫度均為600 K，計算時采用全反射邊界條件。采用蒙特卡羅程序統(tǒng)計特征值和圖8中1號環(huán)板燃料的微觀吸收截面，NECP-X計算條件與環(huán)板燃料柵元問題一致，采用不同方法計算的1號環(huán)板燃料235U和238U微觀吸收截面及其偏差示于圖9。由圖9可見，基于等效幾何理論的共振計算方法最大偏差發(fā)生在238U第23群，為-3.2%，基于MOC的偽核素子群方法最大偏差也在該能群，為-1.79%，而傳統(tǒng)子群方法最大偏差發(fā)生在235U第25群，為5.57%?；贛OC的偽核素子群方法精度最高，基于等效幾何理論次之，傳統(tǒng)子群方法最差。特征值計算結果列于表3，蒙特卡羅程序統(tǒng)計結果為0.404 48±2.0×10-5，3種方法的特征值計算結果相差不大，而基于等效幾何理論的共振計算方法的計算速度為基于MOC的偽核素子群方法的203.6倍，為傳統(tǒng)子群方法的319.77倍。上述微觀截面及特征值計算結果表明，基于等效幾何理論的共振計算方法在環(huán)板燃料的共振計算中具有良好的計算精度，較基于MOC的偽核素子群方法具有非常明顯的加速效果。

表3 環(huán)板燃料組件問題中不同共振計算方法的計算時間及特征值計算結果

圖8 環(huán)板燃料組件幾何

圖9 環(huán)板燃料組件問題中235U和238U吸收截面及偏差的計算結果

2)蜂窩煤狀燃料組件問題

根據(jù)2.1節(jié)中柵元問題建立蜂窩煤狀燃料組件問題，其組件幾何如圖10所示，所有材料溫度均為600 K，計算時采用全反射邊界條件，NECP-X計算條件與蜂窩煤狀燃料柵元問題一致，除統(tǒng)計組件內(nèi)所有柵元的功率外，還統(tǒng)計了圖10中數(shù)字2處燃料柵元(記為柵元2)的235U和238U的微觀吸收截面。采用不同方法計算的燃料柵元2的235U和238U微觀吸收截面及其偏差示于圖11，各方法的計算時間和特征值計算結果列于表4。由圖11可知，基于等效幾何理論的共振計算方法截面計算結果最大偏差為-1.72%，基于MOC的偽核素子群方法計算最大偏差發(fā)生在238U第25群，為1.69%，傳統(tǒng)子群方法截面計算結果最大偏差發(fā)生在235U第25群，為7.47%。從表4的特征值計算結果可看出，基于等效幾何理論的共振計算方法的特征值偏差為0 pcm，基于MOC的偽核素子群方法的特征值偏差為-35 pcm，傳統(tǒng)子群方法的特征值偏差為-171 pcm，從計算時間來看，基于等效幾何理論的共振計算方法相較于基于MOC的偽核素子群方法加速比為380.05，相較于傳統(tǒng)子群方法加速比為211.47，基于MOC的偽核素子群方法較傳統(tǒng)子群方法時間更長的原因在于，基于MOC的偽核素子群方法進行等效計算時需要多次采用MOC搜索二維等效問題，當燃料柵元增多時，基于MOC的偽核素子群方法的計算時間將會進一步增加。3種方法計算的功率分布基本一致，其中基于等效幾何理論的共振計算方法的計算功率與參考解的偏差如圖12所示，燃料柵元功率在靠近導向管柵元及儀表管柵元處較高，柵元功率最大偏差為-0.128%，發(fā)生在組件邊緣處，可看出，由于蒙特卡羅程序統(tǒng)計的偏差，柵元功率存在不對稱的情況，柵元功率RMS為0.043%。截面計算結果、特征值結果和功率分布結果均表明，基于等效幾何理論的共振計算方法在處理蜂窩煤狀燃料組件問題時精度較高，而傳統(tǒng)子群方法截面計算結果較差，與基于MOC的偽核素子群方法相比，基于等效幾何理論的共振計算方法在保證較高的微觀吸收截面計算精度的情況下，計算時間大幅減少，與傳統(tǒng)子群方法相比，基于等效幾何理論的共振計算方法具有更高的計算精度和計算效率。

圖10 蜂窩煤狀燃料組件幾何

圖11 蜂窩煤狀燃料組件問題中235U和238U吸收截面及偏差的計算結果

表4 蜂窩煤狀燃料組件問題中不同共振計算方法的計算時間及特征值計算結果

圖12 蜂窩煤狀燃料組件問題柵元功率及偏差分布計算結果

2.3 JRR-3全堆芯問題

采用本文提出的基于等效幾何理論的共振計算方法對JRR-3全堆芯問題進行模擬計算，幾何與材料參數(shù)參見文獻[17]，采用二維控制棒全提問題進行模擬，模擬中簡化外部輕水反射層的建模，僅建模到重水反射層。

特征值參考解為1.374 20±3.6×10-5，NECP-X特征值計算結果為1.377 58，特征值偏差為338 pcm，相比上述問題，該問題的各項異性更強。JRR-3組件功率計算結果示于圖13，網(wǎng)格功率及偏差分布計算結果示于圖14，其中的標準燃料組件燃料板長度方向等距劃分20個網(wǎng)格，跟隨燃料組件等距劃分16個網(wǎng)格。

由圖13可見，最大組件功率偏差為-1.182%，組件功率分布呈中間低兩邊高的趨勢。由圖14a可見，網(wǎng)格高功率發(fā)生在組件角點附近，由圖14b可見，雖然有些網(wǎng)格的功率偏差超過5%，但大部分功率偏差均小于±3%，網(wǎng)格功率RMS為1.088%?？紤]到板狀燃料研究堆的強各向異性，特征值與功率計算偏差在可接受范圍內(nèi)。

圖13 JRR-3組件功率及偏差分布計算結果

圖14 JRR-3網(wǎng)格功率及偏差分布計算結果

3 結論

針對全局-局部耦合策略中非棒狀燃料的計算問題，提出了基于等效幾何理論的復雜燃料共振計算方法，并將其應用于NECP-X程序中。該方法通過基于復雜幾何燃料逃脫概率、復雜幾何燃料到周圍結構材料的碰撞概率以及丹可夫修正因子的守恒，將復雜幾何燃料等效為規(guī)則的一維棒狀燃料，在等效一維棒狀燃料上利用偽核素子群計算燃料的有效自屏截面。數(shù)值結果表明：針對環(huán)板燃料柵元、蜂窩煤狀燃料柵元、環(huán)板燃料組件、蜂窩煤狀燃料組件和JRR-3全堆芯問題，該方法具有較高的計算精度與計算效率。