劉仕倡，王 侃，陳義學(xué)

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206；2.清華大學(xué) 工程物理系，北京 100084)

蒙特卡羅方法可根據(jù)原子核的反應(yīng)截面來(lái)隨機(jī)模擬中子與不同核素相互作用時(shí)的反應(yīng)概率。材料溫度將影響靶核熱運(yùn)動(dòng)，從而對(duì)反應(yīng)截面產(chǎn)生影響，這就是所謂的多普勒效應(yīng)。在核反應(yīng)堆的高保真模擬中，燃料和慢化劑有非常精細(xì)的溫度變化和分布。因此，對(duì)于不同溫度下的截面必須采用特殊處理。

傳統(tǒng)方法是預(yù)先產(chǎn)生不同溫度點(diǎn)的截面用于插值?；诰€性插值點(diǎn)數(shù)據(jù)卷積的Sigma1方法是預(yù)處理工具中最常用的方法，如NJOY[1]。然而，當(dāng)溫度點(diǎn)數(shù)較多時(shí)，預(yù)生成的方法會(huì)帶來(lái)內(nèi)存占用方面的問(wèn)題。近年來(lái)，為減少內(nèi)存的占用，在多物理耦合計(jì)算中提出了幾種在線溫度截面處理技術(shù)。MCNP6[2]提出了基于級(jí)數(shù)展開的溫度擬合方法。Becker等[3]提出了利用隨機(jī)抽樣算法計(jì)算有效多普勒展寬截面的方法。為提高隨機(jī)抽樣的效率，Serpent程序還提出了另一種拒絕抽樣算法(TMS方法)[4]。在OpenMC[5]中使用了另一種稱為多極點(diǎn)(multipole)表示的方法。Dean等[6]提出了應(yīng)用Gauss-Hermite積分來(lái)計(jì)算多普勒展寬方程的方法。

然而，現(xiàn)有的方法均有其局限性。例如，MCNP6中的溫度擬合方法依賴于預(yù)先生成的系數(shù)。OpenMC中的多極點(diǎn)表示方法依賴于預(yù)先生成的多極點(diǎn)數(shù)據(jù)。拒絕抽樣方法不能使用徑跡長(zhǎng)度估計(jì)器，這將限制該方法的應(yīng)用。其他方法如Sigma1法和傳統(tǒng)的Gauss-Hermite法存在效率低的問(wèn)題。因此，有必要提出一種新的在線多普勒展寬方法，該方法不需進(jìn)行溫度擬合和多極點(diǎn)數(shù)據(jù)生成等額外的處理，從而提高計(jì)算精度和效率。

本文提出一種新的Gauss-Hermite求積方法，并在蒙特卡羅程序RMC[7]中實(shí)現(xiàn)，以提高傳統(tǒng)Gauss-Hermite求積方法的效率。

1 基于積分的多普勒展寬方法

σ(|v-v′|,0)P(v′,T)dv′

(1)

1.1 Sigma1方法

Sigma1方法是利用點(diǎn)截面的線性化特性進(jìn)行卷積，將式(1)轉(zhuǎn)化為式(2)：

[e-α(v-V)2-e-α(v+V)2]dV

(2)

α=M/2kT

(3)

式中，v和V分別為中子速度和中子相對(duì)于靶核的速度。為簡(jiǎn)化式(2)，定義了如下變量：

x2=αV2

(4)

y2=αv2

(5)

將式(2)轉(zhuǎn)化為：

(6)

再定義：

(7)

式(6)可拆分為兩部分：

式(6)最終可轉(zhuǎn)化為一系列Hn的疊加，該函數(shù)定義如下：

(8)

式中，n為次數(shù)。

在實(shí)際應(yīng)用中，式(8)的求解無(wú)論采用誤差余函數(shù)還是Taylor展開均非常耗時(shí)。

1.2 Gauss-Hermite正交基函數(shù)法

Gauss-Hermite正交基函數(shù)定義為：

(9)

為使式(7)形式類似式(9)的形式，定義變量z=x-y，將式(7)寫為：

(z+y)2e-z2dz

(10)

式中：z為求積組的取值點(diǎn)；n為求積組取值點(diǎn)的數(shù)目，如圖1所示；y代表中子速度，x=z+y代表相對(duì)速度。對(duì)于式(10)，注意到指數(shù)項(xiàng)的存在，z的取值范圍為-4≤z≤4即可。需注意的是，式(10)的積分下界為負(fù)無(wú)窮。因此當(dāng)y>4時(shí)，式(10)可用Gauss-Hermite多項(xiàng)式表示，而y<4時(shí)則不能。因此，在應(yīng)用Gauss-Hermite方法時(shí)，必須考慮適用范圍為y>4。y<4時(shí)將采用高斯-勒讓德求積組或Sigma1算法。在Romano的實(shí)現(xiàn)[8]中也采用了相同的處理方法。在Romano的實(shí)現(xiàn)中，使用了16點(diǎn)正交求積組。

圖1 Gauss-Hermite多項(xiàng)式求積組Fig.1 Diagram of Gauss-Hermite quadrature

然而，退回到高斯-勒讓德求積組或Sigma1算法將降低Gauss-Hermite方法的效率。此外，Gauss-Hermite方法在能量網(wǎng)格中搜索相對(duì)能量x，以及插值相應(yīng)的截面時(shí)效率較低。

2 改進(jìn)Gauss-Hermite方法

針對(duì)傳統(tǒng)Gauss-Hermite方法存在的問(wèn)題，提出了改進(jìn)Gauss-Hermite方法。解決這些問(wèn)題的一些關(guān)鍵技術(shù)如下。

2.1 低能段的修正

眾所周知，低能段截面是相對(duì)比較光滑的，如圖2所示。一般地，低能散射截面變化很小，接近1個(gè)常數(shù)；而低能吸收、裂變截面服從1/v率。因此，由式(1)可知，經(jīng)過(guò)多普勒展寬后，吸收截面保持不變，而散射截面應(yīng)乘以修正項(xiàng)，如式(11)所示。

圖2 900 K時(shí)238U總截面Fig.2 Total cross section of 238U at 900 K

(11)

erf是誤差函數(shù)，并且：

(12)

式中，n為核素。

通過(guò)上述處理，對(duì)于低能光滑段，利用式(11)修正散射截面時(shí)，不使用Gauss-Hermite方法，這可減少Gauss-Hermite方法的使用頻率，代之以計(jì)算成本較低的修正項(xiàng)。從而避免了Gauss-Hermite方法在低能段y<4不適用的問(wèn)題。

低能區(qū)與共振能區(qū)的分界點(diǎn)很重要，必須采用合理的方法將低能區(qū)與共振能區(qū)劃分開。

2.2 共振段的開始和結(jié)束能量

如上所述，低能區(qū)與共振能區(qū)的分界點(diǎn)可稱為共振的起始能量或低能光滑段的結(jié)束點(diǎn)。低能光滑段的結(jié)束點(diǎn)可根據(jù)每個(gè)核素所在的最大溫度下多普勒展寬的作用范圍來(lái)確定，從而得到低能光滑段的結(jié)束點(diǎn)。對(duì)于某個(gè)中子能量點(diǎn)，根據(jù)中子與靶核之間的最大相對(duì)能量和最小相對(duì)能量，可找到上界和下界。然后計(jì)算出上、下界的差值，將低能光滑段結(jié)束點(diǎn)作為共振的起始能量，如圖3所示。

圖3 低能光滑段的結(jié)束點(diǎn)Fig.3 End point of low energy smooth segmen

利用上述算法，圖3中低能光滑段結(jié)束點(diǎn)為1.25 eV。而從圖2可看出在4.41 eV附近有1個(gè)很小的共振峰?？梢姡岢龅纳辖?下界算法可有效地找出低能光滑段結(jié)束點(diǎn)。

共振段結(jié)束點(diǎn)的查找方法與NJOY的處理一致，即可分辨共振區(qū)的上界、閾能反應(yīng)對(duì)應(yīng)的最小能量以及1 MeV三者的最小能量。

2.3 多普勒展寬哈希表

另一個(gè)重要的改進(jìn)是多普勒展寬哈希表的使用。在Gauss-Hermite方法中，尋找相對(duì)能量的能量網(wǎng)格序號(hào)非常耗時(shí)。RMC中采用了哈希表方法，來(lái)加速能量網(wǎng)格查找的速度[9]。例如，當(dāng)前中子能量為E，哈希表先根據(jù)E找到對(duì)應(yīng)的能量分箱i，然后再在能量分箱i中使用二分法查找能量E對(duì)應(yīng)的網(wǎng)格序號(hào)。

當(dāng)哈希表確定了能量分箱i后，可得到分箱i的能量上、下界Eupper和Elower，而根據(jù)Eupper和靶核速度可得到最大相對(duì)能量Emax，根據(jù)Elower和靶核速度可得到最小相對(duì)能量Emin。當(dāng)中子能量在哈希表的分箱i內(nèi)時(shí)，多普勒展寬后的相對(duì)能量將位于多普勒展寬哈希表的分箱i內(nèi)，如圖4所示。

圖4 多普勒展寬哈希表Fig.4 Diagram of Doppler broadening Hash map

3 結(jié)果與分析

采用典型17×17壓水堆組件模型，如圖5所示，組件模型參數(shù)列于表1。對(duì)多普勒展寬的精度和效率進(jìn)行了測(cè)試。在每個(gè)柵元中均有燃料、氣隙(16O)、輕水(1H、16O)。燃料溫度為900 K，其他材料溫度為600 K。蒙特卡羅計(jì)算中，每代使用100 000個(gè)中子，共1 000代，其中300個(gè)非活躍代。選取了采用基于ENDF/B-Ⅶ.0的精確溫度庫(kù)的RMC結(jié)果作為參考。

圖5 典型壓水堆組件模型Fig.5 Typical PWR assembly model

參數(shù)數(shù)值燃料密度,g/cm310.196 燃料富集度,%3燃料芯塊半徑,cm0.409 6燃料包殼外徑,cm0.475 0柵距,cm1.26冷卻劑密度,g/cm31.0

改進(jìn)Gauss-Hermite方法采用300 K的ACE庫(kù)作為基礎(chǔ)數(shù)據(jù)庫(kù)。在輸運(yùn)及燃耗計(jì)算中，均采用了哈希表加速，并對(duì)燃料棒中的中子能譜進(jìn)行統(tǒng)計(jì)。

3.1 輸運(yùn)計(jì)算

對(duì)于輸運(yùn)計(jì)算，燃料中僅有235U、238U、16O。表2對(duì)kinf進(jìn)行了比較?？煽闯鰝鹘y(tǒng)Gauss-Hermite方法和改進(jìn)Gauss-Hermite方法與精確庫(kù)的差別均在1倍標(biāo)準(zhǔn)偏差左右。傳統(tǒng)Gauss-Hermite方法的計(jì)算時(shí)間是參考值的9.3倍，而改進(jìn)Gauss-Hermite方法時(shí)間僅增加了10%。

表2 壓水堆組件算例的kinf和計(jì)算時(shí)間Table 2 kinf and calculation time of PWR assembly case

圖6比較了3種方法的中子能譜，可看出傳統(tǒng)和改進(jìn)Gauss-Hermite方法的能譜均與參考值很接近。大部分相對(duì)誤差在5%以內(nèi)。較大的相對(duì)誤差僅出現(xiàn)在通量本身很小的能量點(diǎn)。

圖6 壓水堆組件算例的能譜相對(duì)誤差Fig.6 Relative error of energy sprectrum for PWR assembly case

3.2 輸運(yùn)-燃耗耦合計(jì)算

第2個(gè)算例為燃耗后(58.065 MW·d/kg(HM))的PWR組件，燃料中共134種核素。此次計(jì)算還增加了燃料中核素的燃耗相關(guān)單群截面統(tǒng)計(jì)，因此相當(dāng)于輸運(yùn)-燃耗耦合計(jì)算進(jìn)行到58.065 MW·d/kg(HM)的一步輸運(yùn)計(jì)算。

從表3可看出，改進(jìn)Gauss-Hermite方法的時(shí)間增加了29%，而kinf僅相差6.5 pcm。而傳統(tǒng)Gauss-Hermite方法的時(shí)間為精確庫(kù)的3.45倍。改進(jìn)Gauss-Hermite方法在帶燃耗相關(guān)單群截面統(tǒng)計(jì)的輸運(yùn)計(jì)算的效率較TMS方法(1.98倍)和OpenMC的Multipole方法(1.59倍)的高，因此更適合應(yīng)用于輸運(yùn)-燃耗耦合計(jì)算。

表3 帶單群截面統(tǒng)計(jì)的燃耗后PWR組件算例的kinf和計(jì)算時(shí)間Table 3 kinf and calculation time of burned PWR assembly case with one group cross section

4 結(jié)論

本文提出了一種改進(jìn)Gauss-Hermite方法，并在蒙特卡羅程序RMC中進(jìn)行了編程實(shí)現(xiàn)。壓水堆燃料組件算例的輸運(yùn)和燃耗計(jì)算結(jié)果表明，改進(jìn)Gauss-Hermite方法精度與傳統(tǒng)Gauss-Hermite方法相當(dāng)，而且可大幅提高計(jì)算效率，特別是輸運(yùn)-燃耗耦合計(jì)算的效率。該方法簡(jiǎn)潔而直接，避免了額外的預(yù)處理過(guò)程，為在線多普勒展寬提供了一種新的有效方法。因此，該方法能適用于各種蒙特卡羅粒子輸運(yùn)程序，具有良好的應(yīng)用前景，為高保真多物理耦合計(jì)算奠定基礎(chǔ)。