朱雁凌,郝琛,張樂瑞

(1.哈爾濱工程大學 核安全與仿真技術重點學科實驗室,黑龍江 哈爾濱 150001;2.清華大學 核能與新能源技術研究院,北京 100084)

近年來,精細化中子物理計算是核領域興起的一個技術研究熱點。哈爾濱工程大學核科學與技術學院開發(fā)了全堆芯一步法高保真中子輸運計算程序(High fidelity NEutron transport program for 3-D nuclear reactor whole core,HNET)。HNET基于多級加速理論,在三維并行、多級CMFD計算框架下,采用堆芯一維軸向先進節(jié)塊計算方法或離散縱標法耦合堆芯二維徑向特征線法的2-D/1-D方法開展全堆芯三維精細化中子輸運計算。其中,共振計算得到的有效共振截面直接影響中子輸運計算的精度,是HNET程序的重要計算模塊。最初HNET版本中采用HELIOS程序的47群多群微觀核截面數(shù)據(jù)庫進行共振計算。應用更普適的核數(shù)據(jù)庫,增加HNET計算程序的適用性,需要拓展新的數(shù)據(jù)庫。但前提需要對于共振計算方法進行優(yōu)化,確保基于新的數(shù)據(jù)庫,可高效、精確地開展共振計算。

傳統(tǒng)的共振計算方法主要包括等價理論方法、超細群方法和子群方法等。超細群法[1]將共振能區(qū)劃分為上萬群,可以準確地描述共振能群截面隨能量的變化,但是計算效率較低。等價理論方法[2]具有較高的精度和效率,但是難以處理復雜幾何問題。子群方法[3-4]具有較高的精度以及良好的幾何適應性,在國際上被廣泛應用,如HELIOS[5]、DeCART[6]、nTRACER[7]等程序,因此HNET采用子群方法開展共振計算。但是,傳統(tǒng)子群方法在計算子群參數(shù)時,容易出現(xiàn)子群參數(shù)為負值的問題,并且沒有考慮多共振核素之間的干涉效應,因此需要進一步處理和優(yōu)化。Peng等[8]采用最優(yōu)化擬合法生成子群參數(shù);Joo等[4]采用本底迭代法計算有效共振截面;Bacha等[9]結合子群方法與共振干涉因子法與來處理局部共振干涉效應。此外,傳統(tǒng)子群方法在開展共振計算時,存在頻繁調(diào)用中子輸運求解器以求解子群固定源方程的現(xiàn)象[10],會嚴重影響計算效率。

HDF5文件格式能夠以靈活、高效的I/O功能支持各種海量、復雜的數(shù)據(jù),具有讀寫速度快、便于用戶快速進行數(shù)據(jù)查找等優(yōu)點[11]。本文基于自主開發(fā)的HDF5格式69群截面數(shù)據(jù)庫,研究并優(yōu)化了子群共振計算方法,實現(xiàn)了精細化中子物理計算程序HNET的高效、精確共振計算,為后續(xù)中子輸運計算提供了精確的多群共振截面。針對子群參數(shù)的計算,本文采用帕德近似法計算子群參數(shù),并通過調(diào)整擬合點對其優(yōu)化,使計算獲得的子群參數(shù)符合物理意義,以確保子群參數(shù)計算的準確性。為提高計算效率,本文采用等效單共振群方法,并對共振核素進行分組,只對每組的代表性核素進行固定源方程求解以提高計算效率。最終,基于上述方法,開發(fā)了HNET共振計算模塊,對無限均勻問題、燃料柵元問題及VERA基準題進行計算,并將結果與蒙特卡羅程序統(tǒng)計結果比較分析,驗證了優(yōu)化后共振計算模塊的準確性。

1 理論模型

1.1 子群參數(shù)計算

子群方法的基本思想是將共振能群從截面的角度進行離散,通過子群截面和子群概率來描述劇烈振動的核反應截面值。子群形式的能群有效共振截面可表示為：

(1)

式中：g表示能群;N表示子群數(shù)目;σx,n和pn分別為反應道x的第n個子群截面和子群概率;φn為第n個子群通量。由式(1)可知,有效共振截面只與子群參數(shù)和子群通量相關。子群計算方法的流程主要為：先計算子群參數(shù),然后采用輸運求解器求解子群固定源方程得到子群通量,最后用子群通量歸并子群截面獲得有效共振截面。

應用子群方法開展共振計算時,忽略裂變源項及上散射項,且不考慮共振散射截面,在中間近似條件下,子群固定源方程[4]可表示為：

Ω·?φg,n(r,Ω)+(Σa,g,n(r)+λΣp,g)·

(2)

式中：λ為中間近似因子;Σa,g,n為子群宏觀吸收截面;Σp,g為宏觀勢散射截面,λΣp,g=∑λiNiσp,g,i,Ni和σp,g,i分別為第i個核素的核子密度及其第g群的微觀勢散射截面。從式(2)可以看出,子群通量反映了幾何結構及材料組成對有效共振截面的影響,且子群固定源方程與多群輸運方程在形式上相同,與高保真中子輸運方法相容性好。由于子群固定源方程的源項與通量無關,可同時對多個子群固定源方程獨立求解,提高計算效率。

在制作與背景截面相關的共振截面表時,假設無限均勻介質中只存在一個共振核素R,并忽略非共振核素的吸收截面[12],則由式(2)可得中子標通量計算公式：

(3)

式中：σb為共振核素R的背景截面,其定義為σb=λΣp/NR。當只有共振核素存在時,背景截面達到最小值λRσp,R,通量“凹陷”的程度較深。當背景截面非常大時,通量大小接近于1,此時共振核素對通量的影響很小。將式(3)代入到式(1)中,可建立中間近似條件下有效共振截面與子群參數(shù)和背景截面的關系：

(4)

(5)

式中：cn與dn為與子群參數(shù)相關的多項式,且dN-1的值為1。結合文獻[10],將數(shù)據(jù)庫中2N-1組背景截面及對應的共振截面代入(5)求解擬合系數(shù)cn與dn,進而得到子群吸收截面σa,n與子群概率pn。其他反應道(總截面、裂變截面等)子群截面的計算表達式為：

(6)

式中：ei為多項式系數(shù),可通過最小二乘法或切比雪夫擬合法對不同背景截面下的其他反應道截面進行擬合求得。

如果共振核素R某一能群截面的波動較小,而選取的子群數(shù)目又過多,會導致擬合多項式的階數(shù)過高,出現(xiàn)子群截面為負的現(xiàn)象。為使子群參數(shù)符合物理含義,在計算得到擬合參數(shù)cn與dn后,將式(5)擬合多項式因式分解(以子群數(shù)目N=3為例)：

(7)

若|x1-y1|≤ε或者|x2-y2|≤ε(x1

(8)

若|x1-y1|≤ε且|x2-y2|≤ε(x1

為驗證子群參數(shù)的準確性,將子群參數(shù)代入到式(4),反推數(shù)據(jù)庫中不同背景截面下的共振截面,并將其與數(shù)據(jù)庫比較,計算偏差可表示為：

(9)

為避免出現(xiàn)子群參數(shù)為負值的情況,并保證子群參數(shù)的準確性,在計算子群參數(shù)時,需要不斷地調(diào)整擬合背景截面點,直到計算偏差δerror小于1%且子群參數(shù)均為正值。為增加帕德近似法中背景截面點的選取方案,HDF5格式多群截面數(shù)據(jù)庫提供30個背景截面點。

1.2 子群固定源方程優(yōu)化

傳統(tǒng)的子群方法需要分別對各個共振核素所有共振能群的每個子群截面所對應的子群固定源方程進行求解。為減少求解子群固定源方程的次數(shù),可采用等效單共振群方法進行優(yōu)化[4]。由式(2)可知,子群固定源方程的源項只與散射源項有關,且中間近似因子λ和勢散射截面σp隨能量變化的幅度很小,因此在求解子群固定源方程時,可將所有共振能群等效為一個能群。此外,子群固定源方程與子群概率無關,而與各個子群截面的值有關。在同一區(qū)域,如果不同共振核素的子群截面相同,即使它們的子群概率不同,它們的子群通量也會相同。因此,可以事先挑選一組子群吸收截面σm,只對這組特定子群截面所對應的子群固定源方程進行求解。特定子群截面的個數(shù)和取值應盡可能地代表實際燃料子群截面的范圍,本文選取的10個特定子群截面分別為：0.1、10、50、100、200、300、500、1 000、2 000和10 000 b。

當某一區(qū)域存在多個共振核素時,各共振核素間的共振峰可能相互重疊,形成共振干涉效應。本文結合HELIOS程序[5]采用Bondarenko迭代方法[9]對共振干涉進行修正。同時,為了進一步提高計算效率,對共振核素進行分組以加速Bondarenko迭代方法。如表1所示,基于共振核素共振峰之間是否重疊,將共振核素分為若干組,并在每組中挑出一個代表性核素。認為同一組的共振核素間存在共振干涉效應,不同組間的共振干涉效應忽略不計。

表1 共振核素組別及代表性核素Table 1 Resonance categories and representative isotopes

基于上述方法,在開展共振計算時,只需針對各組代表性核素特定的子群截面進行子群固定源方程求解,減少了求解子群固定源方程所需的時間。以235U和236U為例,當特定子群吸收截面的數(shù)目取為5時,子群固定源方程的計算次數(shù)為5,而傳統(tǒng)子群方法需要的計算次數(shù)為130(2×13×5)。此時,式(2)中的Σa及λΣp可由無限稀釋條件下共振能群吸收截面對應的共振積分Ig,∞加權得到：

(10)

(11)

式中：下標i表示歸屬于第k組的共振核素;下標r表示第k組的代表性共振核素;σm為特定的子群吸收截面。

由于共振截面表是在均勻系統(tǒng)下制作而成,而實際情況多為非均勻系統(tǒng)。因此在求得特定子群通量φm后,通過逃逸截面Σe將非均勻系統(tǒng)的背景截面與均勻系統(tǒng)等價：

σbm=σb+Σe(σm)/NR

(12)

逃逸截面Σe的定義為：

(13)

在求出與σm相應的逃逸截面Σe后,實際子群截面對應的逃逸截面采用對數(shù)插值獲得,然后通過式(12)計算背景截面,最后將背景截面代入到式(3)計算子群通量。結合Bondarenko迭代方法,共振核素i的有效共振吸收截面[4]可表示為：

(14)

式中：σx,i表示其他共振核素j的吸收截面對當前共振核素i的影響：

(15)

由于σx的大小取決于σa的值,而σa的值又與σx相關,因此需要通過迭代計算獲得最終結果。

2 數(shù)值驗證

本文采用基于ENDF/B-VII.0加工得到的HDF5格式69群截面數(shù)據(jù)庫,其共振能群為第15～27群,共振能量范圍為9 118～4 eV。為驗證共振計算結果的準確性,將蒙特卡洛程序的統(tǒng)計結果作為參考值。在蒙特卡洛程序中設置500代計算,每代投入50 000粒子,前50代為非活躍代數(shù)。在獲得蒙特卡洛程序的統(tǒng)計結果后,采用通量體積權重的方法計算出能群平均截面的參考值。

2.1 無限均勻問題

為驗證子群參數(shù)的準確性,如表2所示設置僅由單個共振核素(235U或238U)和1H所組成的無限均勻問題,不考慮實際幾何結構、材料非均勻性等因素。通過調(diào)整核子密度比例來改變共振核素背景截面的大小,其中1H與235U核子密度的比例選為1 000、500、100;1H與238U核子密度的比例選為10、5、1。

表2 無限均勻問題核素組成Table 2 Compositions for infinite homogeneous cases

當溫度為600 K時,無限均勻問題有效共振吸收截面的相對誤差如圖1所示。由圖1(a)可知,235U除了第27群(9.88～4 eV)外,各吸收截面都具有較高的精度,最大相對誤差不超過0.2%。從圖1(b)可以看出,隨著1H與238U核子密度比例的減小,238U背景截面逐漸減小,各能群有效截面的相對誤差在不同程度上都有所增大,尤其是在case 1～6(此時238U背景截面在20 b左右)。其原因主要是在計算子群參數(shù)時,并沒有選取很小的背景截面作為擬合點,導致當238U的背景截面較小時,子群參數(shù)對截面的擬合精度有所下降。

圖1 無限均勻問題吸收截面的相對誤差Fig.1 Relative error of absorption cross sections for infinite homogeneous cases

圖2對比了溫度為293、600、900和1 100 K時,235U(case 1-3)和238U(case 1-5)有效共振吸收截面的相對誤差。由圖2可知,在所選的溫度范圍內(nèi),235U和238U吸收截面的相對誤差在溫度上沒有明顯的趨勢性變化,溫度變化對子群參數(shù)計算精度的影響比較小?？傮w上,在不同背景截面和不同溫度下,單共振核素無限均勻問題有效共振吸收截面的相對誤差均小于1%,這表明子群參數(shù)的準確性高。

圖2 不同溫度下吸收截面的相對誤差Fig.2 Relative errors of absorption cross sections at different temperatures

2.2 燃料柵元問題

為驗證共振計算模塊處理共振干涉效應的能力,針對只含燃料和慢化劑的單柵元問題進行計算,其中燃料棒半徑為0.409 6 cm,柵元半邊長為0.63 cm。如表3所示,將燃料的組成設置為238U、235U和16O,并根據(jù)富集度設置不同的核子密度。慢化劑由1H和16O組成,它們的核子密度分別為0.022 116 3×1024、0.044 232 6×1024cm-3。燃料和慢化劑的溫度均為600 K。

表3 不同富集度下燃料組成Table 3 Fuel composition at different enrichment

圖3顯示了不同富集度下,燃料棒內(nèi)共振核素棒平均吸收截面的相對誤差。棒平均吸收截面是指采用通量-體積加權方式歸并得到的燃料棒活性區(qū)平均截面。

圖3 燃料棒吸收截面的相對偏差Fig.3 Relative errors of absorption cross sections in fuel rod

從圖3可以看出,不論是235U還是238U,棒平均吸收截面的相對誤差在共振能區(qū)低能段比高能段大,這是因為在低共振能群,235U和238U的共振干涉效應較為強烈。此外,238U的相對誤差整體上要大于235U,這說明共振干涉效應對238U的影響更大。

由圖3可知,235U第25群(27.70～15.97eV)棒平均吸收截面的相對誤差最大。一方面,這是因為235U第25群截面存在若干個可分辨的共振峰,有明顯的共振自屏效應,而且238U第25群截面存在孤立的大共振峰,導致這2個共振核素在第25群的共振峰出現(xiàn)嚴重的相互影響。另一方面,這是由于Bondarenko迭代方法在計算一個共振核素時,假設其他共振核素的截面都不隨能量的變化而變化,在處理共振干涉效應時存在一定的誤差[8]。整體而言,除了235U的第25群,235U和238U各個棒平均吸收截面的相對偏差均在3.5%以內(nèi),能夠保持在一個較低的水平。

表4為不同富集度下燃料柵元問題keff的計算結果,表中keff參考值由蒙特卡洛程序提供。當富集度從1.8%增大到3.1%時,有效增殖因子的計算偏差略有增大,但結果仍在可接受范圍內(nèi)?？傮w上看,本文所實施共振計算方法能夠較好地處理共振干涉效應。

表4 不同富集度下燃料柵元問題keff計算結果Table 4 Calculating results of keff for fuel pin cell at different enrichment

2.3 VERA基準題

本文選擇VERA系列基準題[13]中的單柵元問題和單組件問題對HNET程序的共振計算模塊做進一步驗證。從文獻[13]可獲得基準題的幾何結構及相應的材料組成,其中VERA基準題1/4組件的結構如圖4所示,由燃料棒、導向管和儀表管組成。HNET程序相關計算參數(shù)的設置如下：每個卦限選取64個方位角和3個極角,單柵元問題的特征線間距為0.03 cm,單組件問題的特征線間距為0.05 cm。

圖4 VERA基準題1/4組件結構Fig.4 Geometry of 1/4 lattice for VERA benchmark

單柵元或者單組件各自的4個問題具有相同的幾何結構,它們所不同的是燃料溫度、慢化劑溫度以及核子密度。表5和表6分別為VERA單柵元問題和單組件問題有效增殖因子的計算結果,表中keff參考值由文獻[13]提供。由表5和表6可知,各基準題keff計算結果的最大偏差為109 pcm,總體上計算結果與參考值相符。

表6 VERA單組件問題keff計算結果Table 6 Calculating results of keff for VERA lattice problems

圖5顯示了基準題2B中1/8組件歸一化功率相對誤差絕對值的分布情況,功率分布的參考值由文獻[13]提供。由圖5可知最大相對誤差絕對值為3.424%,且誤差較大的燃料棒主要分布在組件外圍。整體上歸一化功率的平均相對誤差為0.925%,功率分布的誤差在可接受范圍內(nèi)。綜合有效增殖因子和功率分布可知,在典型壓水堆常規(guī)問題上,HNET程序的共振計算模塊具有較高的計算精度。

圖5 基準題2B歸一化功率的誤差Fig.5 Error distribution of normalized pin power of Problem 2B

3 結論

1)針對所研究的單共振核素無限均勻問題,有效共振吸收截面的相對誤差均小于1%,這表明優(yōu)化后的帕德近似法可以保證子群參數(shù)的準確性。

2)Bondarenko迭代方法能夠較好地處理不同共振核素間的共振干涉效應。

3)通過典型基準例題的計算分析,驗證了有效增殖因子和歸一化功率的準確性。這表明優(yōu)化后的子群方法在保證計算效率的前提下,具有較高的計算精度,適用于典型壓水堆常規(guī)問題。