吳文斌， 王志強(qiáng)， 鄭競超， 謝擇懌， 趙晨， 彭星杰

(1.中山大學(xué) 中法核工程與技術(shù)學(xué)院，廣東 珠海 519082； 2.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

特征線方法是一種幾何適應(yīng)能力強(qiáng)、精度高、并行潛力大的中子輸運(yùn)計算方法，它已成為組件程序及高保真數(shù)值堆所采用的主流方法，目前已經(jīng)有一批針對柵元呈正方形排列的特征線輸運(yùn)計算程序[1]。然而，針對六角形幾何堆芯布置，相對來說開發(fā)滯后，目前，DeCART程序[2]、SONG輸運(yùn)計算程序[3]、NECP-X程序[4]和ANT-MOC程序[5]等具備針對六角形幾何的特征線輸運(yùn)求解的能力。六角形燃料組件因其排列緊湊，可以最大化冷卻劑流速及其傳熱效率，在先進(jìn)多用途反應(yīng)堆的堆芯設(shè)計中廣泛使用[6]。因此，開發(fā)具備六角形幾何處理能力的中子輸運(yùn)求解器具有重要的工程實(shí)用價值。

OpenMOC程序[7]是MIT開發(fā)和維護(hù)的MOC開源程序，可用于方形堆芯和組件的特征線輸運(yùn)計算。本文在OpenMOC程序的基礎(chǔ)上，開發(fā)了六角形堆芯中子輸運(yùn)計算程序OpenMOC-HEX,實(shí)現(xiàn)六角形模塊化射線布置和基于區(qū)域分解的消息傳遞接口(message passing interface，MPI)并行計算。并采用多個基準(zhǔn)例題對程序進(jìn)行數(shù)值驗(yàn)證和并行效率分析。

1 六角形特征線理論模型

1.1 六角形模塊化射線布置

模塊化射線布置中，特征線僅在小模塊(柵元或組件)內(nèi)生成；小模塊的短特征線通過合理的布置能夠精確連接，構(gòu)成長特征線。模塊化射線布置方法在矩形幾何的MOC程序中廣泛使用[8-11]。本文參考文獻(xiàn)[2]將模塊化射線布置方法應(yīng)用到六角形幾何中。模塊化射線布置必須滿足2個條件：1)模塊內(nèi)短特征線在模塊間精確連接；2)模塊反射邊界上射線與反射線精確連接。而，射線角度、射線間距的選取只有滿足特定的條件，才能使上述2個條件成立。本文射線布置方式為：

2)在(0,π/6)角度范圍內(nèi)方位角數(shù)目N/12，全部方位角分為N/12組；

3)(0,π/6)角度范圍內(nèi)，用戶預(yù)期的方位角平均分布為：

(1)

(2)

式中：a為六角形棱邊長度；ceil為向上取整。

圖1 模塊化射線布置示意Fig.1 Schematic diagram of modular ray arrangement

5)模塊化修正后的每組方位角的第1個角度φs與射線間距δs為：

(3)

6)每組內(nèi)6個方位角的射線間距相等，每組內(nèi)的6個方位角δs為：

(4)

其中，極角求積組選取Tabuchi-Yamamoto求積組[12]。

1.2 六角形柵格幾何處理

為了靈活快速定義六角形組件，需要建立六角形柵格(Lattice)。如圖2所示，每個柵元單元可填充不同的燃料棒。在特征線方法中，需要獲取特征線段所在網(wǎng)格的材料參數(shù)。因此，若給定任意點(diǎn)在笛卡爾坐標(biāo)系下的坐標(biāo)(x,y)，需求出在該點(diǎn)所屬柵元單元的編號，進(jìn)而得到具體的柵元材料信息(燃料棒信息)。六角形柵格幾何處理中，編號求解方法如圖2所示。

圖2 六角形柵格幾何處理Fig.2 Hexagonal lattice geometry processing

(5)

根據(jù)計算得到(i1,i2)編號坐標(biāo)，坐標(biāo)點(diǎn)可能位于(i1,i2)、(i1+1,i2)、(i1,i2+1)、(i1+1,i2+1)4個六角形柵元單元中，需要進(jìn)行進(jìn)一步的判斷。具體判斷方法為：計算坐標(biāo)點(diǎn)(x,y)與4個六角形柵元單元中心點(diǎn)的距離，距離最短者則為坐標(biāo)點(diǎn)(x,y)所在的柵元單元。

1.3 區(qū)域分解并行實(shí)現(xiàn)

如圖3所示，每個MPI進(jìn)程都負(fù)責(zé)處理一個六角形組件，不同組件間通過角通量通信將計算耦合在一起，實(shí)現(xiàn)MPI并行。程序使用非阻塞模式通信，以減少通信的等待時間，避免通信死鎖。為了能將角通量數(shù)據(jù)準(zhǔn)確發(fā)送給目標(biāo)進(jìn)程，需找出相鄰組件的對應(yīng)關(guān)系。

圖3 4圈六角形組件構(gòu)成的堆芯MPI進(jìn)程編號Fig.3 MPI ranks diagram for a core composed of four-ring hexagonal assemblies

本文對于y方向的六角形組件的進(jìn)程編號(rank)從最上方六角形組件開始，沿著順時針方向，從外向里遞增。例如，4圈六角形y方向組件構(gòu)成的堆芯中，MPI進(jìn)程編號(rank)如圖3所示。

不同圈數(shù)下，六角形的鄰邊對應(yīng)關(guān)系不相同。為了找出鄰邊六角形組件的進(jìn)程編號(rank)，先根據(jù)進(jìn)程編號(rank)對應(yīng)的(i1,i2)坐標(biāo)，利用(i1,i2)坐標(biāo)找到對應(yīng)邊的鄰居六角形(i1,i2)坐標(biāo)。如圖4所示，再根據(jù)鄰邊對應(yīng)六角形組件的(i1,i2)坐標(biāo)找出對應(yīng)的進(jìn)程編號(rank)。

圖4 (i1,i2)坐標(biāo)法尋找鄰邊區(qū)域Fig.4 (i1, i2) coordinate method to find neighboring regions

2 特征線程序數(shù)值驗(yàn)證

2.1 C5G7 六角形堆芯基準(zhǔn)題

使用六角形C5G7 2D基準(zhǔn)題對程序進(jìn)行驗(yàn)證，如圖5所示，該堆芯由MOX燃料組件，UO2燃料組件和最外圈的反射組件組成，外邊界為真空邊界條件。

圖5 六角形C5G7 2D基準(zhǔn)題示意Fig.5 Diagram of the hexagonal C5G7 2D benchmark

圖6是UO2燃料組件和MOX燃料組件的具體布置。每個六角形組件邊長為11.9 cm，均采用如圖7所示的柵元規(guī)格，由邊長0.78 cm的六角形組成，柵元包含外徑為0.54 cm的棒，用于填充燃料、裂變腔和導(dǎo)向管。反射組件全部由水組成，組件邊長同為11.9 cm。

圖6 C5G7 2D基準(zhǔn)題燃料組件Fig.6 Hexagonal C5G7 2D benchmark fuel assemblies

圖7 燃料、裂變腔和導(dǎo)向管柵元示意Fig.7 Schematic diagram of fuel, fission chamber and guide tube cell

特征線建模參數(shù)為60個均勻分布方位角，特征線徑向間距為0.03 cm；對于燃料柵元和慢化劑柵元，網(wǎng)格劃分參數(shù)均為徑向2圈、周向4扇區(qū)；求解采用線性源近似；特征值收斂精度1×10-5，裂變率收斂精度10-5。測試平臺采用AMD EPYC 7T83 64-Core Processor 處理器，主頻最大2.45 GHz，內(nèi)存容量512 GB，雙路CPU共128計算核心256線程。

計算值keff為1.162 32，參考值為1.162 44，參考解來自蒙特卡羅程序McCARD[2]，其中功率最大棒的相對誤差為-1.02%；功率最小棒的相對誤差為0.17%；相對誤差絕對值的最大值為1.56；相對誤差絕對值的平均值為0.30；相對誤差均方根為0.36；功率加權(quán)相對誤差為0.26。

六角形C5G7 2D基準(zhǔn)題棒功率分布相對誤差如圖8所示，其中誤差為

(6)

式中P為棒功率值。

圖8 六角形C5G7 2D基準(zhǔn)題棒功率分布相對誤差Fig.8 Relative error of power distribution of hexagonal C5G7 2D benchmark

由上述驗(yàn)證結(jié)果可知，keff計算結(jié)果與參考解相差12×10-5，燃料棒功率分布最大相對偏差為1.56%(出現(xiàn)在堆芯外部功率較低處)，表明程序計算精度良好。

2.2 多圈六角形組件并行效率基準(zhǔn)題

該基準(zhǔn)題由水、UO2和MOX燃料組成，材料的宏觀截面由C5G7基準(zhǔn)題提供。2圈組件堆芯結(jié)構(gòu)如圖9所示，堆芯由2圈共7個相同的組件構(gòu)成。每個組件的中心和外層燃料由大柵元(圖10(a))組成，材料填充為MOX-4.3%，中間層燃料由小柵元(圖10(b))組成，材料填充為UO2，剩余部分用水填充。

圖9 2圈六角形組件并行效率基準(zhǔn)題示意Fig.9 Schematic diagram of parallel efficiency benchmarking of two-turn hexagonal components

圖10 燃料柵元幾何尺寸Fig.10 Fuel geometry

使用該基準(zhǔn)題組件進(jìn)行弱可擴(kuò)展性分析。逐漸增加六角形的組件圈數(shù)進(jìn)行測試，原本為2圈7個組件，增加一圈后為3圈19個組件，其中組件數(shù)Z與MPI進(jìn)程數(shù)一致，對應(yīng)關(guān)系為：

Z=(r-1)×r×3+1

(7)

式中r為圈數(shù)。

依次遞增圈數(shù)進(jìn)行弱可擴(kuò)展性分析測試。測試不啟用OpenMP并行，測試平臺CPU為Intel Xeon Platinum 9242 CPU@2.30 GHz (單節(jié)點(diǎn)共96核)，使用ICC 17.0.5編譯器和Intel 2017 MPI并行庫。測試均為100次迭代。每種圈數(shù)計算4次，計算時間取后3次結(jié)果平均值。不同的圈數(shù)的測試結(jié)果如圖11所示。

圖11 弱可擴(kuò)展性分析運(yùn)行時間隨MPI進(jìn)程數(shù)變化Fig.11 Weak scalability analysis of runtime variation with the number of MPI processes

理想情況下，組件和MPI進(jìn)程數(shù)增加，計算時間應(yīng)不變。實(shí)際情況中，7～91進(jìn)程在單節(jié)點(diǎn)內(nèi)計算，隨著進(jìn)程數(shù)的增加計算時間較明顯地增加；可能原因是單節(jié)點(diǎn)內(nèi)使用核數(shù)的增加使得每個核平均可使用的緩存更少甚至主頻降低，故將91進(jìn)程的運(yùn)行時間作為并行效率計算基準(zhǔn)更為合理。91到1 027進(jìn)程，計算時間只增加約2 s，并行效率為87%(91進(jìn)程計算時間/1 027進(jìn)程計算時間)，表明程序具有良好的并行擴(kuò)展性和并行效率。

3 結(jié)論

1)本文基于開源特征線程序OpenMOC實(shí)現(xiàn)了六角形模塊化射線布置和計算功能，并基于區(qū)域分解策略實(shí)現(xiàn)MPI多進(jìn)程并行，開發(fā)了OpenMOC-HEX程序，使程序能對六角形堆芯進(jìn)行計算。對基準(zhǔn)題的測試表明：程序具有良好的計算精度，六角形C5G7 2D基準(zhǔn)題計算解與參考解keff與參考解相差12×10-5，燃料棒功率分布最大偏差為1.56%；程序具有良好的MPI并行效率，滿足計算精度和擴(kuò)展性上的要求。

2)OpenMOC-HEX程序尚不具備CMFD加速功能，對于C5G7 六角形堆芯基準(zhǔn)題外迭代次數(shù)約500次，計算代價較大。

程序暫未實(shí)現(xiàn)射線軸向的布置，無法計算三維六角形問題。下一步將實(shí)現(xiàn)基于空間區(qū)域分解的六角形CMFD加速功能和三維六角形全堆芯計算功能。