基于EFEN-SP3方法的pin-by-pin中子動力學(xué)計算方法研究

謝偉華，曹良志，李云召

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

pin-by-pin的壓水堆三維瞬態(tài)計算可直接給出安全分析所需的單棒功率密度分布而不需要引入組件均勻化近似和精細功率重構(gòu)帶來的誤差[1]。但pin-by-pin中子動力學(xué)計算以柵元為最大計算網(wǎng)格，一方面使中子擴散計算的假設(shè)不再滿足，另一方面計算量和存儲量也都增加了幾個量級。本文在西安交通大學(xué)自主研發(fā)的壓水堆堆芯燃料管理計算程序系統(tǒng)NECP-Bamboo2.0[2]的基礎(chǔ)上，采用全隱式向后差分方法離散SP3中子動力學(xué)方程組的時間項，采用指數(shù)函數(shù)展開節(jié)塊-SP3(EFEN-SP3)方法[3]求解含裂變介質(zhì)的固定源方程。其中，由于EFEN中采用的平源近似假設(shè)與緩發(fā)中子先驅(qū)核濃度計算時認為的緩發(fā)中子先驅(qū)核濃度形狀和中子通量密度形狀相同不兼容，造成計算誤差會隨時間步累積。因此，本文采用高階源展開的指數(shù)函數(shù)展開方法[4]抑制計算偏差隨時間步不斷增加的趨勢，并采用一系列加速技術(shù)對EFEN-SP3的固定源問題求解過程進行加速。

1 理論模型

三維多群時空SP3中子動力學(xué)方程組經(jīng)全隱式向后差分離散后可得tn+1時刻的固定源方程組[5]：

(1)

其中：

(2)

相應(yīng)的緩發(fā)中子先驅(qū)核濃度方程為：

(3)

若Pn為第n階勒讓德多項式、hi為節(jié)塊在i方向上的長度(cm)，將源項在x、y、z3個方向上用勒讓德多項式展開：

(4)

其中：sn為勒讓德多項式展開系數(shù)；N為多項式展開最大階數(shù)。

基于EFEN-SP3方法，可得到中子通量密度解析解[4]為：

ψ(x,y,z)=Cu+eκx+Cu-e-κx+Cv+eκy+

Cv-e-κy+Cw+eκz+Cw-e-κz+ψ0+

(5)

將其近似為多項式形式，可得：

(6)

其中：Cu±、Cv±、Cw±分別為x、y、z3個坐標系上的展開系數(shù)；κ為自定義系數(shù)，由總截面和擴散系數(shù)求得。

本文采用3種加速方法對EFEN-SP3的固定源問題求解過程進行加速，具體包括Ks因子[6]、LW外推[7]和粗網(wǎng)再平衡[8]。

2 數(shù)值驗證與分析

本文基于建立的pin-by-pin中子動力學(xué)計算方法，研發(fā)了相應(yīng)的計算程序EFEN-K，并對高階源展開方法和加速方法進行了驗證與分析。

2.1 高階源展開方法驗證

二維TWIGL擴散基準題是一簡化的兩群瞬態(tài)問題，堆芯是由3個燃料區(qū)組成的邊長為160 cm的正方形，組件尺寸為8 cm×8 cm，詳細的堆芯幾何布置示于圖1，材料的截面參數(shù)參考文獻[9]。瞬態(tài)過程由燃料區(qū)1的熱群宏觀吸收截面Σa，2發(fā)生如下線性變化引起：

(7)

圖1 二維TWIGL基準題幾何布置Fig.1 Geometric distribution of 2D TWIGL benchmark

設(shè)堆芯的初始歸一化功率為1.0，分別以10-2、10-3、10-4和10-5s的時間步長模擬計算0.5 s的瞬態(tài)過程。

平源近似下的功率曲線如圖2a所示，參考Bamboo-Transient[10]的相對偏差曲線如圖2b所示。由圖2可看出：計算結(jié)果不隨時間步長的縮小而收斂，在時間步長取10-5s時，相對偏差開始顯著增大。高階源展開下的功率曲線如圖3a所示，參考Bamboo-Transient的相對偏差曲線如圖3b所示。由圖3可看出：計算結(jié)果隨時間步長的縮小而逐漸收斂，并穩(wěn)定在0.35%左右，主要是SP3方法與擴散近似之間的方法差異帶來的偏差。

圖2 平源近似下的歸一化功率和功率相對偏差Fig.2 Normalization power and relative deviation of power with flat-source approximation

圖3 高階源展開下的歸一化功率和功率相對偏差Fig.3 Normalization power and relative deviation of power with high-order source expansion

2.2 加速方法驗證

參考VERA基準題[11]設(shè)計壓水堆三維單組件pin-by-pin算例，其徑向右下1/8幾何布置如圖4所示。整個組件包含1個裂變室、24個導(dǎo)向管與264根燃料棒，柵元間距為1.26 cm，最外層?xùn)旁?.04 cm的水隙，軸向半組件高度為1 250 cm，上部有20 cm的反射層，之外是真空邊界條件，下部采用全反射邊界條件。4群截面參數(shù)列于表1，動力學(xué)參數(shù)列于表2。瞬態(tài)過程由導(dǎo)向管材料的總截面在初始時刻發(fā)生如下階躍變化引起，近似模擬0.02 s內(nèi)控制棒插入過程。

圖4 加速方法驗證算例徑向1/8示意圖Fig.4 Sketch of 1/8 verification case for acceleration method

(8)

采用4種不同方案進行計算：1) 不加速；2) Ks因子加速；3) Ks因子和LW外推(Ks+LW)加速；4) Ks因子、LW外推和粗網(wǎng)再平衡(Ks+LW+CM)加速。細網(wǎng)計算的網(wǎng)格劃分是徑向每個柵元1個計算網(wǎng)格，軸向網(wǎng)格為5.0 cm。粗網(wǎng)再平衡的網(wǎng)格合并方案為：徑向x方向6個網(wǎng)格和y方向6個網(wǎng)格與z方向5個網(wǎng)格合并成1個網(wǎng)格。

表1 加速方法驗證算例單組件截面參數(shù)Table 1 Cross section of fuel rod in verification case for acceleration method

續(xù)表1

表2 加速方法驗證算例動力學(xué)參數(shù)Table 2 Kinetics parameter of verification case for acceleration method

注：β=2.133 33×10-4pcm，λ=0.08 s-1

針對4種方案的計算結(jié)果比較如圖5a所示，以無加速下的結(jié)果為參考解，另外3種方案結(jié)果的相對偏差如圖5b所示。由圖5可見，加速方法對結(jié)果影響很小，可忽略。對4種方案的計算總時間進行統(tǒng)計，結(jié)果列于表3，可看出：Ks加速可取得47.6倍左右的加速；LW外推在此基礎(chǔ)上可獲得1.99倍的加速；粗網(wǎng)再平衡可在此基礎(chǔ)上獲得1.42倍的加速；3種加速技術(shù)共可實現(xiàn)134.9倍的加速。

圖5 4種方案SP3計算的歸一化功率與相對偏差Fig.5 Normalization power and relative deviation of power of SP3 calculation with 4 cases

方案計算時間/s加速比無加速99 675.388Ks加速2 092.70747.6Ks+LW加速1 050.05894.9Ks+LW+CM加速739.055134.9

2.3 三維多組件問題

該算例參照C5G7-TD基準題[12]，徑向幾何如圖6所示，由兩個UOX組件與兩個MOX組件相隔布置而成，四周為全反射邊界，活性區(qū)高度為42.68 cm，頂部為21.34 cm的軸向反射層，外為真空邊界條件，底部無反射層，為全反射邊界條件。柵元均勻化少群常數(shù)由壓水堆柵格計算程序Bamboo-Lattice2.0提供。瞬態(tài)過程中，所有材料的總截面與散射截面的線性變化，用于模擬硼濃度的變化。具體為在0～1 s內(nèi)所有能群的總截面線性增加0.1%，所有能群的散射截面線性減少0.1%；在1～2 s內(nèi)截面線性恢復(fù)初始截面。瞬態(tài)模擬過程為3 s。

圖6 pin-by-pin問題徑向幾何布置Fig.6 Radial configuration of pin-by-pin problem

EFEN-K采用SP3近似，計算時徑向網(wǎng)格為1個柵元(1.26 cm×1.26 cm)1個網(wǎng)格，軸向網(wǎng)格大小為5.355 cm。參考解為采用輸運動力學(xué)程序NECP-X[13]進行同樣pin-by-pin計算得到的結(jié)果。

歸一化功率隨時間的變化如圖7所示，可看出：SP3近似的結(jié)果與參考解存在一定偏差，最大偏差在1.0 s時刻，歸一化功率分別為0.138與0.124，主要是因為歸一化功率絕對值非常低導(dǎo)致相對偏差較大。另外，參考解采用高階輸運MOC方法與SP3近似之間存在誤差，在實際計算時SP3近似方法需結(jié)合合適的柵元均勻化技術(shù)。取第0.25、0.75、1.875、2.0和3.0 s時刻的徑向二維功率分布和軸向一維功率分布，以NECP-X的功率分布結(jié)果為參考解，得到SP3近似的徑向棒功率均方根相對偏差和最大相對偏差以及軸向一維功率分布的最大相對偏差和均方根相對偏差，結(jié)果列于表4?？煽闯觯簭较蜃畲笙鄬ζ詈途礁鄬ζ罹霈F(xiàn)在0.75 s，分別為-2.598%與1.006 4%；軸向最大相對偏差和均方根相對偏差均出現(xiàn)在第0.75 s，分別為0.389 9%和0.121 6%。其中，第0.75 s的徑向棒功率相對偏差分布如圖8所示，可看出相對偏差呈對角布置，MOX組件多為負偏差，UOX組件多為正偏差，但最大相對偏差均不超過3%。取第0.75 s的軸向功率分布及其相對偏差列于表5，可看出相對偏差主要在靠近反射層的活性區(qū)位置，但均不超過0.5%。

圖7 三維多組件問題歸一化功率隨時間的變化Fig.7 Normalization power change with time of 3D multi-assembly problem

時刻/s徑向最大相對偏差/%徑向均方根相對偏差/%軸向最大相對偏差/%軸向均方根相對偏差/%0.25-0.995 60.368 00.338 10.113 90.75-2.5981.006 40.389 90.121 61.8750.734 40.275 30.338 60.114 32.0-0.895 40.291 80.321 20.111 83.00.901 50.291 00.301 90.108 8

圖8 三維多組件問題第0.75 s的SP3近似徑向棒功率相對偏差分布Fig.8 Relative deviation distribution of radial rod power of SP3 calculation for 3D multi-assembly problem at 0.75 s

表5 三維多組件問題第0.75 s軸向棒功率分布及其相對偏差Table 5 Axial rod power distribution and relative deviation of 3D multi-assembly problem at 0.75 s

3 結(jié)論

為實現(xiàn)壓水堆三維pin-by-pin中子動力學(xué)計算，本文研究了基于全隱式向后差分方法和EFEN-SP3方法的中子動力學(xué)計算技術(shù)，通過高階源展開方法消除了固定源問題求解過程中的平源近似與緩發(fā)中子先驅(qū)核濃度計算中的高階近似之間的不兼容問題，通過一系列加速方法有效提高了計算效率，并通過一系列算例的計算進行了數(shù)值驗證與分析。數(shù)值結(jié)果表明：1) pin-by-pin中子動力學(xué)計算能直接給出安全分析所需的單棒功率密度分布，且計算結(jié)果與中子輸運時空動力學(xué)程序NECP-X的pin-by-pin計算結(jié)果符合良好；2) 高階源展開技術(shù)可有效抑制計算偏差隨時間步的累加效應(yīng)；3) 加速技術(shù)可將SP3動力學(xué)計算的求解速度提高134.9倍。