張宏博，趙 晨，彭星杰，趙文博，陳 長(zhǎng)，李 慶，于穎銳，宮兆虎，曾 未，劉 琨，饒俊杰，王 博

(中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

近年來(lái)，借助于高性能計(jì)算平臺(tái)的長(zhǎng)足進(jìn)步，數(shù)字化反應(yīng)堆的高保真、高分辨率建模與多物理模擬計(jì)算能力，對(duì)核能系統(tǒng)的設(shè)計(jì)改進(jìn)與機(jī)理認(rèn)知顯現(xiàn)出重要且特殊的意義。其中，高保真中子學(xué)數(shù)值模擬技術(shù)作為數(shù)字化反應(yīng)堆體系的核心組成之一，成為業(yè)內(nèi)研究熱點(diǎn)。近年來(lái)，僅在以特征線方法(MOC)為基礎(chǔ)的非均勻輸運(yùn)理論體系框架下，就陸續(xù)涌現(xiàn)出了DeCART[1]、nTRACER[2]、MPACT[3]、NECP-X[4]和STREAM[5]等代表性軟件。

SHARK程序是由中國(guó)核動(dòng)力研究設(shè)計(jì)院正在研發(fā)的基于全堆芯確定論非均勻輸運(yùn)理論體系的高保真中子學(xué)計(jì)算軟件。該軟件采用構(gòu)造實(shí)體幾何(CSG)的建模方式，從多群數(shù)據(jù)庫(kù)的截面與共振數(shù)據(jù)出發(fā)，采用改進(jìn)子群方法刻畫有效共振截面的復(fù)雜非均勻效應(yīng)[6]，采用一步化的二維/一維(2D/1D)或準(zhǔn)三維(Quasi-3D)MOC開展堆芯非均勻輸運(yùn)計(jì)算[7]。SHARK程序采用面向?qū)ο蟮脑O(shè)計(jì)理念，采用Python與C++混編開發(fā)，具有良好的可維護(hù)性和用戶友好性。

本文介紹SHARK程序的研發(fā)情況，并采用一系列基準(zhǔn)題對(duì)該程序進(jìn)行測(cè)試驗(yàn)證。

1 方法模型

1.1 多群數(shù)據(jù)庫(kù)

SHARK程序現(xiàn)階段采用基于ENDF/B-Ⅵ.8核評(píng)價(jià)庫(kù)的45群中子能群結(jié)構(gòu)多群數(shù)據(jù)庫(kù)[8]，這種能群結(jié)構(gòu)既能較為充分地捕捉壓水堆堆芯內(nèi)的中子能譜特點(diǎn)，又不會(huì)為堆芯規(guī)模的高保真計(jì)算帶來(lái)過(guò)重的計(jì)算負(fù)擔(dān)，符合程序的應(yīng)用堆型定位。SHARK程序數(shù)據(jù)庫(kù)包含各類核素共計(jì)280余種，基本覆蓋了程序目標(biāo)應(yīng)用對(duì)象的燃料、慢化劑/冷卻劑、毒物、結(jié)構(gòu)材料及相關(guān)燃耗鏈核素，考慮了錒系核素、裂變產(chǎn)物核素、可燃吸收體核素的裂變、俘獲、(n，2n)/(n，3n)和衰變等主要反應(yīng)類型。

1.2 全堆芯共振算法

在子群方法中，有效共振截面σx的表達(dá)形式為：

(1)

式中：wx，n和σx，n為x類型子群參數(shù)；下標(biāo)a為吸收類型，n為子群編號(hào)；φn為子群通量。

在式(1)中，若要求取σx，必須求解具有輸運(yùn)方程形式的子群方程以獲得子群通量。為提高堆芯層面共振計(jì)算效率，引入慢化方程近似解將式(1)改寫為：

(2)

(3)

(4)

式中：σb，n為子群背景截面；λ為中間共振因子；Σp為宏觀勢(shì)散射截面；Σe，n為宏觀等效截面；N為共振核素核子密度。

利用式(2)～(4)計(jì)算有效共振截面的關(guān)鍵是獲得Σe，n，而直接求解同樣是不經(jīng)濟(jì)的。由于Σe與吸收截面σa呈現(xiàn)相對(duì)平緩的單調(diào)關(guān)系，故可構(gòu)造Σe(σa)函數(shù)關(guān)系線性對(duì)數(shù)插值表，以減少子群方程求解頻次。

子群等效截面插值表Σe(σa)是與共振能群和共振核素相關(guān)的。根據(jù)核素共振峰分布及目標(biāo)堆芯的特點(diǎn)，可采用共振核素分類及共振能區(qū)分段兩種策略削減計(jì)算成本。

基于上述簡(jiǎn)化，以MOC求解形如式(5)的子群方程，獲得子群通量并完成制表：

λΣp+(1-λ)Σpφm(r)m=1,…,M

(5)

式中：Ω為方向單位向量；r為中子向徑；ψm為子群m的角通量。

當(dāng)同一材料區(qū)域內(nèi)含有多種共振核素時(shí)，采用迭代思想處理多核素共振干涉效應(yīng)。

對(duì)于共振散射截面，可根據(jù)邦達(dá)連科迭代收斂后的有效共振吸收積分，利用二分法反向搜索背景截面，并據(jù)此計(jì)算有效共振散射截面。

對(duì)于非共振核素，采用移出截面修正因子方法相對(duì)準(zhǔn)確刻畫輕核的慢化效應(yīng)，以考慮共振所致的能譜變化引發(fā)的包殼、慢化劑等材料的散射截面變化。

1.3 非均勻輸運(yùn)算法

SHARK程序包含2D/1D MOC和Quasi-3D MOC兩種非均勻輸運(yùn)求解器。總體思想上，二者都是將三維輸運(yùn)問(wèn)題轉(zhuǎn)化為軸向若干層相互耦合的二維MOC問(wèn)題進(jìn)行求解的。目前程序主要采用2D/1D MOC。

在2D/1D MOC中，三維中子輸運(yùn)方程經(jīng)軸向和徑向積分，可得基本耦合方程為：

(6)

(7)

SHARK程序采用空間有限差分格式的離散縱標(biāo)方法(SN)求解軸向一維方程(式(7))，以充分考慮軸向的非均勻輸運(yùn)效應(yīng)。

對(duì)于2D/1D MOC數(shù)值計(jì)算的潛在不穩(wěn)定性問(wèn)題，采用各向同性泄漏近似和泄漏項(xiàng)分割方法進(jìn)行彌補(bǔ)。

而在Quasi-3D MOC中，通過(guò)引入軸向差分格式近似，將三維中子輸運(yùn)方程轉(zhuǎn)化為二維耦合方程：

(8)

(9)

由式(8)、(9)可見(jiàn)，Quasi-3D MOC不含有橫向泄漏項(xiàng)，方程右項(xiàng)恒正，因此較2D/1D MOC具有更好的數(shù)值穩(wěn)定性[7]。

1.4 加速與并行計(jì)算技術(shù)

對(duì)于大規(guī)模計(jì)算，加速與并行計(jì)算技術(shù)是十分必要的。加速技術(shù)方面，SHARK程序采用最優(yōu)擴(kuò)散粗網(wǎng)有限差分(odCMFD)加速方法[9]，該方法較好地改善了傳統(tǒng)粗網(wǎng)有限差分(CMFD)方法的數(shù)值穩(wěn)定性。

并行策略方面，SHARK程序采用基于空間區(qū)域分解的MPI+OpenMP并行方案，即采用MPI技術(shù)處理空間區(qū)域分解并行，各子區(qū)分別獨(dú)立進(jìn)行特征線追蹤計(jì)算，并與內(nèi)邊界相鄰子區(qū)通信(圖1)；采用OpenMP技術(shù)對(duì)特征線產(chǎn)生和追蹤過(guò)程中的部分循環(huán)體進(jìn)行并行。這種混合并行策略可有效分散計(jì)算與存儲(chǔ)負(fù)荷，具有良好的并行效率。

圖1 SHARK程序空間區(qū)域分解并行策略

綜合以上方法模型，圖2給出程序的計(jì)算流程和基本邏輯關(guān)系。輸運(yùn)計(jì)算最終以特征值和平源區(qū)裂變率為收斂判定準(zhǔn)則。

圖2 SHARK程序基本計(jì)算流程

2 數(shù)值結(jié)果

SHARK程序目前已具備定態(tài)微觀問(wèn)題的中子學(xué)計(jì)算能力，各項(xiàng)測(cè)試驗(yàn)證均在有序推進(jìn)。本文給出了部分基于國(guó)際著名基準(zhǔn)題的計(jì)算結(jié)果。

2.1 宏觀BEAVRS堆芯問(wèn)題

該問(wèn)題改編于裝載有193盒UO2組件的大型商用壓水堆堆芯[7]。堆芯中布置有1.6%、2.4%、3.1% 3種富集度的燃料組件。堆芯活性區(qū)高為380 cm，并在軸向兩端各有40 cm的鐵水反射層(圖3)。

圖3 BEAVRS堆芯幾何示意圖

該問(wèn)題參考解來(lái)自蒙特卡羅程序RMC，以驗(yàn)證SHARK程序輸運(yùn)求解器。其中，SHARK程序采用的活性區(qū)軸向?qū)痈邽?.5 cm；在每個(gè)象限，采用8個(gè)方位角和3極角Yamamoto最佳求積組；特征線寬度為0.02 cm；特征值和源項(xiàng)的收斂準(zhǔn)則分別為1 pcm和10-4。BEAVRS堆芯問(wèn)題計(jì)算結(jié)果列于表1，棒功率分布偏差如圖4所示。由于高保真計(jì)算節(jié)點(diǎn)資源限制，以136核并行為基準(zhǔn)，對(duì)比204、408及578個(gè)CPU核心數(shù)量，計(jì)算時(shí)間及并行效率列于表2。由計(jì)算結(jié)果可見(jiàn)，對(duì)于普通商用壓水堆堆芯，非均勻輸運(yùn)具有較為良好的精度及效率。

表1 BEAVRS堆芯問(wèn)題計(jì)算結(jié)果

表2 BEAVRS堆芯問(wèn)題并行計(jì)算效率

圖4 BEAVRS堆芯問(wèn)題徑向棒功率相對(duì)偏差分布

2.2 二維組件問(wèn)題(VERA#2)

VERA#2為二維單組件基準(zhǔn)題。該基準(zhǔn)題較為全面地考慮了燃料富集度、溫度、可燃毒物、控制棒等常見(jiàn)因素給商用壓水堆組件帶來(lái)的非均勻效應(yīng)。

表3列出VERA#2基準(zhǔn)題的計(jì)算統(tǒng)計(jì)結(jié)果，其中參考解來(lái)自蒙特卡羅程序KENO-Ⅵ[10]。SHARK程序在每個(gè)象限采用16個(gè)方位角和3極角Yamamoto最佳求積組；特征線寬度為0.03 cm(含IFBA問(wèn)題特征線寬度為0.005 cm)；特征值與源項(xiàng)的收斂準(zhǔn)則分別為1 pcm和10-5。對(duì)于特征值結(jié)果，SHARK程序的最大偏差為-220 pcm，大多數(shù)算例的偏差均在±150 pcm以內(nèi)。對(duì)于棒功率分布結(jié)果，絕大多數(shù)算例的最大相對(duì)偏差(MAX)不足±0.30%，均方根偏差(RMS)不足0.10%；對(duì)于含有強(qiáng)吸收體的算例(典型為控制棒)，最大偏差水平有所增大。由此可見(jiàn)，對(duì)于二維單組件系列問(wèn)題，程序的計(jì)算精度是合理和穩(wěn)定的。

表3 VERA#2基準(zhǔn)題計(jì)算結(jié)果

2.3 三維組件問(wèn)題(VERA#3)

VERA#3基準(zhǔn)題是一個(gè)壓水堆三維單組件問(wèn)題，包含無(wú)可燃毒物(3A@600 K)和含有16根Pyrex毒物(3B@565 K)兩個(gè)算例。這兩個(gè)算例改編自真實(shí)的商用壓水堆組件。除徑向布置，該基準(zhǔn)題力求還原燃料組件在軸向上的端塞、壓緊彈簧、管座、底座和定位格架等結(jié)構(gòu)細(xì)節(jié)。因此該問(wèn)題可較為充分地反映軟件對(duì)壓水堆軸向非均勻效應(yīng)的模擬能力。

同樣以KENO-Ⅵ為參考，并默認(rèn)采用相同的軸向?qū)觿澐?，?列出SHARK程序計(jì)算VERA#3基準(zhǔn)題的宏觀統(tǒng)計(jì)結(jié)果，圖5示出VERA#3基準(zhǔn)題徑向棒功率偏差分布，圖6示出VERA#3基準(zhǔn)題軸向功率分布。由表4可見(jiàn)，特征值的計(jì)算偏差較小，徑向棒功率偏差與二維組件問(wèn)題的偏差水平相一致(圖5)。

算例：a——3A；b——3B

表4 VERA#3基準(zhǔn)題計(jì)算結(jié)果

對(duì)于軸向功率分布，SHARK程序的最大相對(duì)偏差僅有2.61%，且出現(xiàn)在頂部靠近反射層的低功率區(qū)域。同時(shí)，軸向功率平均偏差(AVG)和均方根偏差也說(shuō)明了軸向功率分布的模擬結(jié)果是準(zhǔn)確的。由圖6可見(jiàn)，軸向功率形狀與KENO-Ⅵ程序契合，且格架效應(yīng)得到了準(zhǔn)確體現(xiàn)。

算例：a——3A；b——3B

2.4 三維多組件問(wèn)題(VERA#4)

VERA#4基準(zhǔn)題是一個(gè)壓水堆組件3×3排布問(wèn)題，在徑向和軸向上都有著比前述算例更強(qiáng)的非均勻效應(yīng)。該系列問(wèn)題包含3個(gè)二維算例和12個(gè)三維算例。本文對(duì)這些算例進(jìn)行了精細(xì)建模與驗(yàn)證計(jì)算。

1)二維算例

根據(jù)中心位置RCCA組件插入情況不同，二維算例包括控制棒全提(ARO)算例4A、銀銦鎘(AIC)控制棒算例4B和碳化硼(B4C)控制棒算例4C。

表5列出VERA#4二維多組件基準(zhǔn)題計(jì)算結(jié)果。由表5可見(jiàn)，SHARK程序能較準(zhǔn)確地捕捉徑向多組件的非均勻效應(yīng)，特征值偏差不超過(guò)-43 pcm。圖7示出棒功率相對(duì)偏差分布，少數(shù)臨近控制棒位置的棒功率分布相對(duì)偏差最大達(dá)到2%左右(相對(duì)功率低至0.500以下)，絕大多數(shù)棒功率相對(duì)偏差不超過(guò)±1%。

表5 VERA#4二維多組件基準(zhǔn)題計(jì)算結(jié)果

圖7 VERA#4二維多組件基準(zhǔn)題徑向功率相對(duì)偏差分布

2)三維算例

與VERA#3基準(zhǔn)題相同，VERA#4基準(zhǔn)題同樣考慮了軸向的精細(xì)結(jié)構(gòu)。此外，每根控制棒在軸向上被分為兩段，底部為AIC吸收體，上部則為B4C吸收體。該問(wèn)題共包含12個(gè)分別處于不同控制棒位的算例，側(cè)重于比較特征值、控制棒微分價(jià)值(DRW)、控制棒積分價(jià)值(IRW)和軸向功率偏移(AO)。這些參數(shù)能較為全面地刻畫程序?qū)β史植嫉挠?jì)算精度水平。

表6列出VERA#4三維多組件基準(zhǔn)題的宏觀物理參數(shù)計(jì)算結(jié)果。由表6可見(jiàn)，特征值偏差不超過(guò)100 pcm。對(duì)于棒價(jià)值，各棒位處DRW與IRW結(jié)果均吻合良好(圖8)。同時(shí)，AO數(shù)據(jù)與參考解也十分契合。

表6 VERA#4三維多組件基準(zhǔn)題計(jì)算結(jié)果

圖8 VERA#4三維多組件基準(zhǔn)題控制棒價(jià)值

3 總結(jié)

本文介紹了中國(guó)核動(dòng)力研究設(shè)計(jì)院數(shù)字化反應(yīng)堆高保真中子學(xué)程序SHARK的基本研發(fā)情況。SHARK程序采用匹配目標(biāo)堆芯物理特性的多群數(shù)據(jù)庫(kù)、共振計(jì)算與輸運(yùn)計(jì)算方法模型，同時(shí)應(yīng)用了CSG建模、MPI+OpenMP并行等先進(jìn)技術(shù)手段。數(shù)值結(jié)果顯示，SHARK程序已具備定態(tài)微觀中子學(xué)問(wèn)題計(jì)算能力，對(duì)于商用壓水堆相關(guān)基準(zhǔn)題具有良好的計(jì)算精度和效率。根據(jù)項(xiàng)目規(guī)劃，后續(xù)將對(duì)SHARK程序持續(xù)開發(fā)完善，陸續(xù)完成精細(xì)燃耗、多物理耦合、時(shí)空動(dòng)力學(xué)等方法模型與代碼的開發(fā)，并逐步推進(jìn)更加深入和廣泛的測(cè)試驗(yàn)證工作。