“華龍一號”堆芯中子學計算的蒙特卡羅源收斂問題

申鵬飛,霍小東,邵 增,楊海峰,張 鵬,王 侃

(1．中國核電工程有限公司,北京,100840;2.清華大學 工程物理系,北京,100084)

蒙特卡羅算法常用于堆芯高精度建模和高保真物理計算中。蒙特卡羅臨界計算通常采用源迭代過程,使用用戶自定義的源作為初始源,進行輸運模擬并將其后代存入中子庫中;下一代抽樣時,從上一代的中子庫中進行隨機抽樣的方式得到新的源分布,進行輸運模擬并重新獲得中子庫;以此往復,得到收斂源分布后再進行各物理量的統計。

蒙特卡羅的源迭代過程存在源收斂問題,表現為源慢收斂、方差低估計等問題[1-3]。其中,源慢收斂問題體現為從用戶自定義的源分布逐步迭代消除初始源分布的高階項,最終得到本征源分布的過程,在高占優比的模型中耗時較久,甚至需數百非活躍代。而壓水堆堆芯模型是典型的高占優比模型,占優比可達0.98～0.99。方差低估計問題體現為在蒙特卡羅活躍代計算中,統計量的計算方差小于真實方差的現象。造成方差低估計現象的原因是蒙特卡羅每個活躍代的統計結果都被視為相互獨立。在HOOGENBOOM全堆模型中,組件功率的方差低估計系數(under estimation ratio, UER)可大于90%[4]。綜上,在大尺度的高占優比系統中,源收斂問題會顯著影響計算資源的分配和計算結果的分析解讀,因此需針對各個模型進行蒙特卡羅源收斂問題研究。

“華龍一號”(HPR1000)是我國自主研發的第三代核電技術堆型[5],反應堆采用 177 個燃料組件構成的堆芯(簡稱“177 堆芯”),具有完全自主知識產權[6],這一新的堆芯設計方案會導致不同的源收斂特性。為深入分析“華龍一號”堆芯蒙特卡羅臨界計算的源收斂特性,本文使用keff、香農熵及WD(Wasserstein distance)等判據分析了源慢收斂問題,給出了非活躍代數的建議值,同時計算了柵元、組件和全堆3個層次的集總量和分布量的UER。使用組件功率的實測值與計算值進行了對比,并使用半定量預估消除(semi-quantitative variance underestimation elimination,SeVUE))算法計算了代間相關性影響長度,給出了解決源收斂問題的相關計算參數。計算平臺使用國產堆用蒙特卡羅軟件RMC。

1 堆芯建模及源收斂診斷

HPR1000堆芯由 177 個先進燃料組件構成。堆芯活性段高度(冷態)為365.8 cm,堆芯高徑比為1.13。首循環堆芯燃料組件按235U的富集度分3區裝載,富集度分別為1.80%,2.40%,3.10%[6]?！叭A龍一號”首循環堆芯簡化模型如圖1所示。

(a)Top view

(b)Side view

為檢驗源收斂診斷特性,設置不同位置、不同每代粒子數的初始源,并使用keff、香農熵與WD判據進行收斂性分析[7-8]。初始源位置如圖1(b)所示,初始源的參數如表1所列。

表1 初始源的參數

圖2為“華龍一號”堆芯keff隨代數的變化關系。由圖2可見,數十代內,源實現收斂。香農熵是常用的蒙特卡羅源收斂判據,“華龍一號”堆芯香農熵隨代數的變化關系如圖3所示。由圖3可見,不同初始源設置下收斂特性有所不同：當初始源設置在中部時收斂最快,約100代完成收斂;當初始源設置在頂部、底部時收斂較慢,約200代完成收斂;不同每代粒子數對收斂特性沒有影響。

圖2 “華龍一號”堆芯keff隨代數的變化關系

圖3 “華龍一號”堆芯香農熵隨代數的變化關系

Guo等[8]提出的WD判據是根據裂變源分布的數據特征進行收斂診斷,是一種更保守的收斂判據。該判據分別對3個方向的粒子坐標進行排序并計算其WD(xWD,yWD,zWD),不同初始源設置下堆芯WD隨代數的變化關系如圖4所示。由圖4可見,與香農熵類似,不同初始源設置下收斂特性有所不同：當初始源設置在中部時收斂最快,約200代完成收斂;當初始源設置在頂部、底部時收斂較慢,約300～350代完成收斂。此判據下,還觀察到每代粒子數不同時,源收斂特性的差異,當每代粒子數增大時,所需收斂代數略有增加。這是因為每代粒子數增大時,WD距離在收斂后的波動減小,使WD距離收斂到穩態的代數增加。對5種方案的源收斂特性進行綜合比較發現,非活躍代參數應設置為大于300代,且keff判據和香農熵判據均存在提前判斷收斂的現象。

(a)Scheme index 1(top)

(b)Scheme index 2(middle)

(c)Scheme index 3(bottom)

(d)Scheme index 4(top)

(e)Scheme index 5(top)

2 數值計算方差低估計現象

計算了集總量keff、中子代時間、總功率及組件與柵元功率分布的方差低估計現象。蒙特卡羅臨界計算參數如表2所列。參考解及真實方差由40組獨立計算得到,其中每次計算設置不同的隨機數種子,舍棄單次計算方差σa,使用統計方法得到40組計算值的真實方差σr。方差低估計系數ηUER可定義為

(1)

表2 蒙特卡羅臨界計算參數

集總量的方差及方差低估計系數如表3所列。由表3可知,3個參數的方差低估計現象均不明顯,方差低估計系數小于±10%以內,與參考解方差的固有波動相近。

表3 集總量的方差及UER

在“華龍一號”的臨界分析中,組件及柵元功率等分布量更值得關注?！叭A龍一號”堆芯組件功率分布的參考解及方差低估計現象如圖5所示。由圖5(d)可見,組件功率分布的UER為77%～93%,堆芯內部組件的UER相對較低,外側組件較高。原因是堆芯內部中子自由程有限,外側組件本代源分布與上代源分布的相關性高于中心區域組件。全部組件的平均UER為89.1%。

(a)Reference power

(b)Reference relative deviation

(c)Relative deviation of Monte Carlo single calculation

(d)UER

選取了C3,E5,H8,K11,M13等斜對角分布的5個組件(如圖1所示)分析柵元功率分布的方差低估計現象。其中,E5組件的柵元功率分布參考解及方差低估計現象如圖6所示。

(a)Reference power

(b)UER

由圖6(b)可見,柵元功率分布的UER無明顯規律,為23%～53%。5個組件柵元功率分布的平均UER如表4所列。與組件功率分布的方差低估計現象類似,堆芯組件的柵元功率分布的方差低估計現象較不顯著,外圍組件的柵元功率分布的平均方差低估計系數較高,可達40%以上。

表4 5個組件柵元功率分布的平均UER

因此,功率分布等分布量的方差低估計現象需引起足夠關注。在“華龍一號”組件功率分布的計算處理中,應將計算標準差乘以10;在不同位置組件的柵元功率分布計算處理中,通常應將計算標準差乘以1.5。

3 組件實測功率分布對比

將“華龍一號”熱態滿功率(hot full power, HFP)狀態下組件功率的實測值進行比較分析。比較C3,E5,H8組件的實測功率與蒙特卡羅中子學計算的功率差異,并實際觀察蒙特卡羅方差低估計現象。各組件蒙特卡羅計算相對偏差如表5所列。

表5 各組件蒙特卡羅計算的相對偏差

C3,E5,H8組件的實測功率與蒙特卡羅40組單次計算功率分布如圖7所示。由圖7可見,在40組蒙特卡羅獨立計算中,單次蒙特卡羅計算得到的組件功率大多數不能在3倍標準偏差內與實測結果吻合。修正了方差低估計現象后,蒙特卡羅計算的參考解與實測功率在3倍真實標準偏差內吻合良好。因此,在實際堆芯的中子學計算中,須采取適當的措施修正方差低估計現象。

(a)C3

(b)E5

(c)H8

4 代間相關性長度計算

使用SeVUE方法計算“華龍一號”堆芯的代間相關性長度,并消除各個層面的低估計系數[9]。SeVUE方法首先計算蒙特卡羅裂變源分布的平均SWD(Sliced Wasserstein distance)隨相隔代數的變化關系,判斷出該模型的代間相關性長度,并從理論上證明了模型的代間相關性長度就是組統計方法中的最佳組長度。其中,組統計方法是一種將蒙特卡羅活躍代歸并以減弱蒙特卡羅方差低估計的方法,須設置組長度l,將活躍代每隔l代歸為多個組,以組為單位統計物理量并計算標準偏差。由于組間裂變源分布的相關性明顯弱于代間,因此組統計方法可有效減弱蒙特卡羅方差低估計的現象。但組長度l的設置是一個難點,如l設置過大,組數偏少,統計結果的標準偏差的波動就較大;如l設置過小,組間裂變源分布的相關性就可能較大,從而無法完全消除方差低估計現象。

對于“華龍一號”堆芯模型,其裂變源分布平均SWD隨相隔代數的變化關系如圖8所示。由圖8可見,代間相關性長度為裂變源分布平均SWD達到最大值并開始波動的間隔長度,應為約500代。

圖8 “華龍一號”模型裂變源分布平均SWD隨相隔代數的變化關系

SeVUE方法的第二步為蒙特卡羅組統計計算。設置組長度為500代,其他參數與表1相同,得到組件功率分布的UER,如圖9所示。由圖9可見,使用SeVUE方法消除方差低估計現象后,組件功率分布的UER無明顯規律,在-8.9%～14.4%之間波動,平均值為6.3%。

圖9 使用SeVUE方法后組件功率的UER分布

柵元功率分布的方差低估計現象被消除。5個組件的柵元平均UER如表6所列。以E5組件為例,其柵元功率分布的UER如圖10所示。

表6 使用SeVUE方法后,5個組件柵元功率分布的平均UER

圖10 使用SeVUE方法后,E5組件柵元功率分布的UER

由圖10可見,柵元功率分布的UER也無明顯規律,在-25%～23%之間波動,平均值為6.8%。

因此,除對分布量的計算結果進行系數處理外,可使用組長度為500代的組統計計算,可完全消除臨界計算方差低估計現象。

5 結論

本文使用keff分析、香農熵診斷及WD診斷算法分析了“華龍一號”堆芯的蒙特卡羅源慢收斂問題,給出收斂所需的建議非活躍代數為30代。研究了蒙特卡羅方差低估計問題,說明了keff、中子代時間及總功率等集總量的方差低估計現象不顯著;對于功率密度分布,進行了柵元和組件兩個層次功率分布的UER計算,平均UER分別約為30%和90%,且堆芯外側的組件功率和柵元功率方差低估計現象較顯著。與實測值對比,蒙特卡羅單次計算的組件功率大多數不能在3倍標準差內與實測值吻合;修正了方差低估計后,蒙特卡羅計算的參考解可在3倍真實標準差內與實測值吻合。

本文還使用SeVUE算法計算了“華龍一號”堆芯模型的代間相關性長度,確定約為500代。建議使用系數修正或組長度為500代的組統計法進行“華龍一號”堆芯模型的臨界計算分析。