基于子群方法的雙重非均勻性共振計算方法研究

李 頌，張 乾,*，張志儉，趙 強，梁 亮，梁越超，張晉超，婁 磊，李滿倉

(1.哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001；2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610213)

隨著新型壓水堆燃料的發展，雙重非均勻性逐漸成為了堆芯計算的難點。傳統燃料只包含燃料棒與慢化劑之間的一重非均勻性，而全陶瓷微密封(FCM)燃料中的TRISO以顆粒形式彌散分布在碳化硅基體中[1]，組成了第二重非均勻性。由于彌散顆粒具有隨機排布的性質，傳統共振方法無法進行雙重非均勻性系統的共振計算。

為解決此問題，國際上主要提出了3種方法。第1種方法為丹可夫方法[2]，通過中子逃脫概率計算丹可夫因子來修正有效共振截面，主要用于高溫氣冷堆，不適用于FCM燃料球的多層結構。第2種方法為缺陷因子方法[3-6]，通過顆粒與基體的等效均勻化，將問題簡化為傳統一重非均勻系統。該方法簡單高效，但無法給出顆粒內部的精細通量與截面分布。第3種方法為Sanchez-Pomraning方法[7]，可直接求解雙重非均勻性系統特征線方法(MOC)輸運方程，給子群共振計算提供了可行性。然而，子群方法難以處理雙重非均勻性系統中的共振干涉，傳統Bondarenko方法精度較低，干涉因子方法需建立復雜的雙重非均勻性積分表，偽核素方法則需求解大量超細群方程，計算效率較低。為精確計算雙重非均勻性系統的共振截面，本文提出耦合Sanchez-Pomraning方法的改進的子群方法(ISSP)。

1 理論模型

1.1 基于Sanchez-Pomraning方法的MOC模型

傳統MOC無法直接對FCM燃料隨機分布的彌散顆粒進行特征線掃描。為解決此問題，Sanchez-Pomraning方法提出了“更新方程”的概念。首先，通過迭代法計算基體和顆粒混合的等效總截面，如式(1)所示。基于等效總截面，平源區的出射通量如式(2)所示， 其中φas為漸進通量，計算公式如式(3)所示。平源區等效源項q=Σtφas。

(1)

(1-e-ΣtL)(φas-φin)

(2)

(3)

式(2)建立了傳統MOC與雙重非均勻性系統間的聯系。對不含彌散顆粒的平源區，φas計算方法與傳統MOC相同。在彌散燃料區域，基體通量φm可通過直接求解式(1)～(3)得到。基于基體通量，顆粒內部各層通量則通過式(4)求得。

(4)

1.2 改進的子群方法

與傳統子群方法[8]相比，改進的子群方法采用精細能群結構來規避共振干涉處理。子群參數包括子群截面σi及子群權重ωi，可通過中間共振近似及帕德近似方法計算得到[9]。有效共振截面σx與子群參數間的關系為：

(5)

其中：g為當前能群；φi為子群通量，需求解子群固定源方程(式(6))得到。

(6)

其中：Ω和r分別為角度和位置變量；σt′,i=σt,i-(1-λ)σs,i，σt，i和σs，i分別為子群總截面和子群散射截面；Nr為共振核素r的核子密度；λn為核素n的中間共振因子；Σp，n為勢散射截面。式(6)右側為源項，左側第2項為移出項。

在雙重非均勻性條件下，子群移出項、源項與通量計算需要式(1)～(3)來進行修正。為減少固定源方程求解次數，本文結合HELIOS[10]和nTRACER[11]程序對子群固定源方程的優化，提出了等效單群固定源優化方法[9]。將所有共振能群等效為1個共振群，子群固定源方程根據等效共振群的子群截面進行計算，并通過插值得到實際子群通量。等效能群的子群截面計算方法由文獻[9]提供，本文選取取值范圍為10-24～10-20cm2的8個等效子群截面，可滿足不同類型問題計算。

為保證輸運計算效率，本文采用中子慢化方程來獲得共振能群能譜。其移出項由細群總截面計算得到，源項來自上游能群的下散射。求出第1共振群的通量后，下游共振群的源項即可根據散射矩陣和上游能群通量逐群求得。第1共振群的散射源項來自于快群，其能譜由漸進通量1/E近似提供。在雙重非均勻性條件下，慢化方程移出項、源項與通量計算同樣需要式(1)～(3)來修正，并通過式(7)歸并得到最終有效多群截面。

(7)

1.3 基于Sanchez-Pomraning方法的超細群方法

本文選取的超細群在共振能量段共有34 000群，中子與燃料核素碰撞后所能獲得的最大對數能降遠大于超細群的勒寬，因此在每個超細群能量段內中子最多只會碰撞1次。定義fk,n為中子飛越k個超細群后與核素n碰撞的概率，如式(8)所示，每個超細群的源項可由式(9)計算得到。

(8)

(9)

其中：N為核素總數；K為中子與核素n碰撞后能飛越的最大能群數；Δug為第g群勒寬；u0、u、u′為能量變量；Σn，s為散射截面；α=(A-1)2/(A+1)2，A為相對原子質量。

由式(8)可得fk=e-Δuk-1fk-1，將此關系代入式(9)可得超細群源項，如式(10)所示。因此，在求出第1共振群的通量后，下游能群源項即可逐群求解固定源方程獲得。第1共振群的源項由漸進通量近似求解，結合式(1)～(3)求得所有超細群通量后，再根據給定寬群能群結構通過式(7)進行并群計算。

(10)

2 數值驗證

本文使用ENDF-Ⅶ版數據庫，采取HELIOS-1.11程序的47群能群結構[10]，其共振能量范圍為1.855～9 118 eV。本文多群數據庫由開源程序NJOY-2016[13]制作得到，問題基準解由開源蒙特卡羅程序OpenMC[14]提供。蒙特卡羅建模采用顯式隨機模型，粒子總數設置為100 000個，粒子代數為1 600代，前300代不計入統計。

2.1 FCM單柵元問題

1) 不同填充率

圖1 FCM燃料幾何結構示意圖Fig.1 Scheme of geometry structure of FCM fuel

單柵元FCM燃料幾何結構如圖1所示。本算例采用富集度為17.8%的碳化鈾燃料，各材料及幾何信息由文獻[15]提供，系統溫度為300 K。本文依次選取1%、5%、10%、20%、30%和40%的填充率進行計算，有效增殖因數計算結果列于表1，并與基于Sanchez-Pomraning方法的超細群方法(UFGSP)的計算結果進行對比。

從表1可見，ISSP與UFGSP計算結果與基準解一致，各填充率偏差均小于100 pcm，而ISSP的效率更高。同時，表1分別列出了ISSP子群固定源和慢化方程兩個模塊的計算時間以及共振和本征值計算的總時間。子群時間在各填充率下均低于總時間的4%，而慢化方程時間最高占比為28.2%，因此ISSP的并群效率仍有可提升空間。以填充率10%和40%為例，圖2示出不同填充率下238U和235U吸收截面的計算精度。

表1 不同填充率下keff的計算偏差Table 1 Calculating deviation of keff for different packing fractions

圖2 不同填充率下的238U(a)和235U(b)吸收截面的相對偏差Fig.2 Relative deviation of absorption cross-section of 238U (a) and 235U (b) under different packing fractions

ISSP對40%填充率的238U吸收截面的相對偏差均在1%以內，對10%填充率的最大相對偏差上升至1.9%。由于偏差較大能群的截面遠小于其他共振群，因此對最終結果影響很小。UFGSP同樣有上述偏差變化趨勢，且所有能群相對偏差小于1%。整體而言，ISSP和UFGSP均可較高精度地提供各填充率下的共振截面，而ISSP具有更高的計算效率。

2) 毒物問題

本文采用如圖3所示的兩種毒物顆粒，第1種是在TRISO燃料顆粒內增加1層毒物材料(QUADRISO)，第2種是毒物以單獨毒物球(BISO)的形式與TRISO顆粒一起彌散在基體中。本文采用B4C和Gd2O3兩種毒物材料，其中單顆粒問題的QUADRISO填充率為30%，雙顆粒問題中的TRISO顆粒填充率為33.8%，BISO為10.8%。系統溫度為300 K，材料核子密度與幾何尺寸來自文獻[15-16]。

a——QUADRISO；b——BISO圖3 可燃毒物顆粒的幾何結構Fig.3 Geometry structure of burnable poison

B4C毒物在QUADRISO和BISO兩種形式下的238U吸收截面的相對偏差如圖4所示，ISSP和UFGSP均具有較高的計算精度，所有能群的相對偏差均小于1%。QUADRISO的相對偏差偏大，主要是因為與燃料內核毗鄰的毒物層導致顆粒內部的通量波動更為劇烈。含Gd毒物中的155Gd和157Gd對反應率影響較大，其吸收截面的相對偏差如圖5所示。

圖4 含B4C毒物燃料QUADRISO和BISO顆粒238U吸收截面的相對偏差Fig.4 Relative deviation of 238U absorption cross-section in QUADRISO and BISO with B4C poison

155Gd共振峰主要集中在10 eV以下能量段，該范圍能群內ISSP和UFGSP對QUADRISO和BISO的最大截面相對偏差均低于2%。在29.02～47.81 eV及130.1～2 034.7 eV之間，155Gd的最大相對偏差增至3.3%左右，但由于此能量段截面很小，因此對全局反應率的影響也較小。ISSP和UFGSP對157Gd的計算精度和155Gd類似，其最大相對偏差在2.6%左右。毒物問題的keff計算精度列于表2，其偏差均在150 pcm以內。總體來看，ISSP和UFGSP均可對毒物材料進行準確的共振計算。

圖5 155Gd(a)和157Gd(b)吸收截面的相對偏差Fig.5 Relative deviation of absorption cross-section of 155Gd (a) and 157Gd (b)

表2 毒物問題keff的偏差Table 2 Deviation of keff for burnable poison problem

3) 燃耗問題

以填充率為30%的算例為初始條件計算FCM燃料的燃耗過程，其中每個燃耗步的核子密度由OpenMC提供，并由ISSP和UFGSP進行共振和輸運計算。有效增殖因數的偏差如圖6所示，ISSP和UFGSP在大多數燃耗步的偏差在100 pcm以內。兩種方法在燃耗過程初始計算精度相近，隨燃耗深度的增加，共振核素大量累積導致干涉效應復雜，計算的偏差會有所增大，在90 GW·d/tU之后偏差逐漸平穩。從整體上看，ISSP和UFGSP均可較好地計算FCM燃料的燃耗過程。

圖6 燃耗過程中keff計算的偏差Fig.6 Deviation of keff during burn-up process

2.2 組件問題

1) 常規FCM燃料組件

本文選取了13×13二維燃料組件，其中包括160根FCM燃料棒和9根導向管，如圖7所示。其中，FCM燃料選自2.1節填充率為30%的算例，導向管材料與幾何組成參考文獻[15]，系統溫度為300 K。以圖7所標燃料棒1和2為例，其238U與235U吸收截面的相對偏差如圖8所示。

圖7 13×13 FCM燃料組件幾何結構Fig.7 Geometry structure of 13×13 FCM fuel lattice

由圖8可見，238U吸收截面的相對偏差在大部分能群小于1%，最大相對偏差分別為-1.54%和-1.42%。235U的計算精度較高，ISSP的最大相對偏差為1.10%，UFGSP在所有能群均小于1%。有效增殖因數和計算時間列于表3，ISSP和UFGSP計算的偏差分別為16 pcm和-1 pcm，二者均可對FCM燃料組件進行精確計算，而ISSP的計算時間遠小于UFGSP。

圖8 13×13組件燃料棒1、2的238U(a)和235U(b)吸收截面相對偏差Fig.8 Relative deviation of absorption cross-section of 238U (a) and 235U (b) in fuel pin 1 and 2 of 13×13 lattice

表3 FCM燃料組件問題keff的偏差Table 3 Deviation of keff for FCM fuel lattice

2) 含毒物FCM燃料組件

本文構造了4×4含毒物FCM燃料組件，其中包括12根常規FCM燃料棒、2根導向管和2根毒物棒，如圖9所示。本算例毒物棒選自2.1節含Gd毒物的QUADRISO燃料，其余材料及幾何均與13×13算例相同，系統溫度為300 K。圖9所標燃料棒1和2的238U以及燃料棒1中155Gd吸收截面的相對偏差如圖10所示。

相比于常規FCM燃料組件，ISSP和UFGSP對毒物組件238U吸收截面計算的相對偏差有所增加，其最大偏差為2.59%，對大多數能群仍保持1%左右，這是由于通量在毒物棒附近劇烈變化所導致。155Gd截面精度也具有相似的規律。該組件本征值計算結果列于表3，相比于常規組件，該組件有效增殖因數偏差有所增加，但仍在可接受范圍內。同時ISSP具有更高的計算效率。

圖9 含Gd毒物4×4 FCM燃料組件幾何結構Fig.9 Geometry structure of FCM 4×4 fuel lattice with Gd poison

圖10 4×4組件燃料棒1、2的238U(a)和棒1的155Gd(b)吸收截面相對偏差Fig.10 Relative deviation of absorption cross-section of 238Uin fuel pin 1 and 2 (a) and 155Gd in fuel pin 1 (b) of 4×4 lattice

3 結論

為解決雙重非均勻性系統的共振計算問題，提出了改進的子群方法。該方法通過精細共振能譜來規避共振干涉處理，并采用等效單群固定源優化來提高計算效率。通過調用基于Sanchez-Pomraning方法的MOC輸運模塊分別計算子群固定源方程和慢化方程，得到FCM燃料顆粒和基體等材料的多群有效截面，最后再調用該輸運模塊進行本征值計算。數值結果表明，ISSP可充分考慮雙重非均勻性的影響，能精確并高效地計算FCM燃料多種工況下的共振截面。