NECP-X中基于改進泄漏項分割技術的二維/一維耦合輸運方法

趙 晨，劉宙宇，曹良志，吳宏春

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

當前計算機發(fā)展水平下，二維/一維耦合輸運因兼顧計算精度和計算效率已被廣泛應用于全堆芯一步法輸運計算[1]。二維/一維耦合輸運方法將一個三維問題分為一系列徑向二維問題和軸向一維問題，并通過泄漏項進行耦合。國際上，對于二維/一維耦合方法進行了較為廣泛的研究，基于該方法開發(fā)了一系列程序，包括CRX-3D[2]、DeCART[3]、nTRACER[4]、MPACT[5-6]等。

國內外對于二維/一維耦合輸運方法的研究主要集中在精度、穩(wěn)定性、效率、內存幾方面。一維輸運方法研究[7]和泄漏項重構方法[8]提高了二維/一維耦合輸運方法的計算精度；通過基于空間、角度、特征線的多重并行很好地解決了內存及效率問題[9]。由泄漏項造成的穩(wěn)定性問題，成為二維/一維耦合輸運方法中最重要的問題之一。對于徑向二維特征線法(MOC)方程，方程右端的軸向泄漏項造成總源項可能出現負值，由此可能造成計算過程中產生負通量導致迭代發(fā)散。對于軸向采用擴散方法、泄漏項采用各向同性近似的二維/一維耦合方法，密歇根大學通過引入松弛因子提高了穩(wěn)定性[5-6]，但該方法對于各向異性泄漏項并不完全適用，且松弛因子方法并未從本質上解決負源項造成的收斂性問題。針對基于各向同性泄漏的二維/一維耦合方法的不穩(wěn)定問題，MPACT通過泄漏項分割方法，提高了穩(wěn)定性的同時保持了較好的精度。但對于各向異性泄漏項的二維/一維耦合方法，該方法會帶來十分嚴重的精度損失。針對各項異性泄漏項的二維/一維耦合輸運方法，國際、國內均未做相關研究。

本文基于西安交通大學開發(fā)的數值反應堆物理計算軟件NECP-X[10-14]，對原有的基于標通量的泄漏項分割方法進行分析，針對各向異性泄漏項的二維/一維耦合方法提出改進的泄漏項分割方法，在內存稍有增加的情況下，解決二維/一維耦合方法中存在的穩(wěn)定性問題，同時克服原有泄漏項分割方法精度損失的問題。

1 理論模型

1.1 二維/一維耦合輸運方法及負源項問題

角度能量離散后的三維中子輸運方程可寫為如下形式：

Qg(x,y,z)

(1)

其中：ε、η、μ分別為x、y、z方向的方向余弦；Σt,g(r)為總截面；ψg,m(x,y,z)為空間(x,y,z)上m方向的角通量，即中子輸運方程中的未知量；Qg(x,y,z)為總源項，包括裂變源及散射源，具體形式如下：

(2)

其中：χg為裂變譜；(νΣf)g′為第g′群的裂變截面；Σs0,g′→g為第g′群到第g群的散射截面；φg′(x,y,z)為標通量。

將式(1)在軸向層L內進行積分，得到徑向二維方程：

(3)

(4)

將式(1)在徑向棒p內進行積分，得到軸向一維方程：

(5)

(6)

NECP-X程序軸向計算采用SN差分方法，基于不作近似的各向異性泄漏項進行耦合，二維/一維迭代直至特征值和裂變率收斂。迭代過程中，由于軸向泄漏項的存在造成式(3)右端源項可能出現負值，在MOC計算過程中，負源項有可能造成沿特征線的出射角通量、平均角通量為負，從而在二維計算中產生負通量，導致迭代發(fā)散。

1.2 泄漏項分割方法

泄漏項分割方法的基本公式如下：

Σt,g(x,y)ψg,m(x,y)+ΣL,g(x,y)ψg,m(x,y)=

(7)

其中，

各向同性泄漏項和基于標通量的泄漏項分割方法被應用于MPACT程序中，并取得了較好的效果。但對于各向異性泄漏項，基于標通量的泄漏項分割方法會造成不可接受的精度損失，因此需要一種適用于各向異性泄漏項的改進的泄漏項分割方法。

1.3 改進的泄漏項分割方法

式(7)中最精確的選擇是選擇各平源區(qū)各角度的角通量進行泄漏項分割，如式(8)～(10)所示。這樣處理后的式(7)與式(3)完全等效，未作任何近似。然而這種方法需存儲每個平源區(qū)每個角度每個能群的角通量，對于大規(guī)模問題這種方法帶來的內存負擔是無法接受的，因此不具備實際應用價值。

(8)

Qg(x,y),0)

(9)

(10)

對此，本文提出了改進的泄漏項分割方法消除二維/一維耦合方法中的負源項問題，并在內存可接受的情況下保持二維/一維耦合計算的精度。因此應基于角通量的泄漏項分割技術，并且不保存平源區(qū)的角通量從而降低內存需求。

為數值上顯示棒內平均角通量和平源區(qū)角通量的分布特點，設計了1個強泄漏算例，如圖1所示。燃料棒在徑向為3×3布置，沿軸向布置2層燃料和1層水，分析的3個平源區(qū)分布也如圖1所示。計算條件為：1個卦限內采用8個幅角3個極角，第1群的棒內平均角通量分布如圖2所示，第1群的3個平源區(qū)內角通量分布如圖3所示。

圖1 強泄漏算例幾何及棒內平源區(qū)分布

圖2 第1群棒內平均角通量分布

如圖2、3所示，每條線表示特定極角方向上角通量沿幅角的分布，棒內均勻化的角通量分布與3個平源區(qū)的形狀類似，因此棒內均勻化的角通量分布可近似表示棒內各平源區(qū)的各向異性效應?；谝陨辖?，各平源區(qū)內的角通量分布可采用標通量和棒內均勻化的角通量分布重構得到，即：

(11)

其中：m為角度序號；i為平源區(qū)序號；pin代表棒內均勻化通量。

圖3 第1群各平源區(qū)角通量分布

基于這種重構方法，僅需存儲棒內均勻化的角通量，存儲量與泄漏項相同，在可接受范圍內。與原有泄漏項分割方法相比，在較小內存需求下獲得平源區(qū)的角通量分布，實現基于角度各向異性的泄漏項分割技術，提高二維/一維耦合方法計算精度。

2 數值結果

2.1 計算精度分析

1) 強泄漏算例

該算例中，強泄漏和強各向異性效應通過設置徑向真空邊界條件引入，幾何如圖1所示。燃料棒在徑向3×3布置，沿軸向布置2層燃料和1層水?；谶@個3×3算例，可將問題拓展至7×7、17×17等更大規(guī)模，建立不同的各向異性問題，從而研究各向異性對計算精度的影響。

該算例中，二維/一維耦合方法未出現負源項引起的迭代發(fā)散問題，因此可采用不帶泄漏項分割方法的二維/一維耦合方法的結果作為基準解，比較傳統(tǒng)基于標通量的泄漏項分割方法、改進的泄漏項分割方法、基于平源區(qū)角通量的泄漏項分割方法的計算結果。特征值計算結果列于表1。傳統(tǒng)基于標通量的泄漏項分割方法嚴重影響計算精度，泄漏越強各向異性越強，計算精度越差。采用改進的泄漏項分割方法可顯著地提高計算精度，最精確的是基于平源區(qū)角通量的泄漏項分割方法。

表1 強泄漏算例的特征值計算結果

2) C5G7基準題全插棒算例

C5G7 OECD/NEA是知名輸運計算基準題，包括不帶棒(unrodded)、半插棒(rodA)、全插棒(rodB)3個算例，表2列出了全插棒算例的計算結果。計算過程中，燃料棒采用48個平源區(qū)、反射層棒采用64個平源區(qū)，1個卦限內采用8個幅角3個極角，特征線寬度為0.03 cm。軸向上，燃料區(qū)分12層，反射層區(qū)分6層，每層3.57 cm。軸向細網為1 cm。

從表2可看出，改進的泄漏項分割方法特征值與基準解相比相差45 pcm，全局最大棒功率偏差為1.94%，計算精度較好，與不用泄漏項分割方法的計算精度相當。然而基于傳統(tǒng)標通量的泄漏項分割方法，特征值偏差達到300 pcm，最大棒功率偏差為15%。結果表明，傳統(tǒng)泄漏項分割方法對于C5G7這類具有較強各向異性的算例會帶來較為嚴重的精度損失，而改進的泄漏項分割方法可獲得較高的計算精度。

表2 C5G7 rodB算例計算結果

3) VERA-3A基準題算例

VERA基準題發(fā)布于2014年，提供了從柵元到堆芯、從二維到三維的一系列基準題，數據多來自Batts Bar壓水堆核電站，因此該基準題是具有真實參數的實際問題。VERA-3A是一個三維組件問題，在NECP-X程序中，對復雜幾何進行了詳盡描述，包括格架、上下管座、下端塞等。燃料棒劃分40個平源區(qū)、反射層劃分64個平源區(qū)，1個卦限包括8個幅角3個極角，特征線線寬為0.03 cm，軸向劃分66層，采用66個核并行。表3列出了特征值和棒功率計算結果，圖4～6示出了軸向、徑向棒功率分布結果。與KENO計算結果相比，特征值偏差小于20 pcm，最大棒功率偏差約1%，表明改進的泄漏項分割方法具有較高的計算精度。

表3 VERA-3A基準題計算結果

圖4 VERA-3A基準題軸向積分功率偏差

圖5 VERA-3A基準題不用泄漏項分割方法徑向積分功率偏差

圖6 VERA-3A基準題改進的泄漏項分割方法徑向積分功率偏差

2.2 內存分析

本文針對2個7群算例分析了改進的泄漏項分割方法的內存要求，1個是7群軸向6層單組件算例，另1個是C5G7基準題全插棒算例。計算條件與C5G7基準題完全一致，內存使用量列于表4。計算結果表明基于平源區(qū)角通量的泄漏項分割方法占用的內存是不用泄漏項分割方法的8倍，對于C5G7基準題，該方法已無法在64G單節(jié)點上進行計算，因此對于更大規(guī)模問題，這種基于平源區(qū)角通量的泄漏項分割方法是不實用的。另一方面，改進的泄漏項分割方法雖然內存占用有所增加，但增加量較小，對于C5G7基準題僅增加9.5%，該方法在更大規(guī)模問題上仍是可接受的。

表4 基于7群算例的內存分析結果

2.3 穩(wěn)定性分析

在1.1節(jié)介紹了二維/一維方法中潛在的負源項造成迭代發(fā)散的問題。對于VERA-3A基準題算例，計算條件與VERA-3A基準題完全相同，使用outflow輸運修正，分別利用24個核和66個核進行并行計算。因為在迭代過程中，反射層區(qū)域出現了很負的源項，導致了迭代發(fā)散的問題。而改進的泄漏項分割方法可解決負源項造成的迭代發(fā)散問題，計算結果列于表5。

表5 對于VERA-3A基準題的穩(wěn)定性分析結果

3 總結

二維/一維耦合輸運方法中由于泄漏項的存在，導致潛在不穩(wěn)定的問題。針對各向異性泄漏項，傳統(tǒng)的基于標通量的泄漏項分割方法在計算強各向異性問題時，會造成較大的精度損失。本文提出了一種改進的泄漏項分割方法，在不造成精度損失的情況下解決了二維/一維耦合輸運方法中的穩(wěn)定性問題。

改進的泄漏項分割方法中，通過棒內均勻化角通量分布和平源區(qū)標通量在線重構平源區(qū)角通量分布。數值結果表明棒內均勻化的角通量分布可近似模擬平源區(qū)的角度各向異性，從而在不保存平源區(qū)角通量的情況下保持二維/一維耦合輸運方法的計算精度。

通過強泄漏算例、C5G7基準題、VERA-3A基準題分析了改進的泄漏項分割方法的計算精度、內存要求和穩(wěn)定性。數值結果表明，改進的泄漏項分割方法可在增加有限內存的條件下，保持二維/一維耦合方法的計算精度，并顯著提高其穩(wěn)定性。