異構(gòu)并行的高階散射特征線方法及其在臨界實驗裝置模擬中的應(yīng)用

鄒 航,陳 瑩,張 乾,曹 巍,張晉超,梁 亮,宋佩濤,劉 杰

(1.國防科技大學(xué) 并行與分布處理重點實驗室,湖南 長沙 410073;2.高端裝備數(shù)字化軟件湖南省重點實驗室,湖南 長沙 410073;3.中國原子能科學(xué)研究院,北京 102413;4.浙江大學(xué) 物理學(xué)系 浙江近代物理中心 先進(jìn)核能理論與應(yīng)用實驗室,浙江 杭州 310030;5.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院 核安全與仿真技術(shù)國防重點實驗室,黑龍江 哈爾濱 150001;6.西安核創(chuàng)能源科技有限公司,陜西 西安 710077;7.中國輻射防護(hù)研究院,山西 太原 030006)

在輕水堆中,各向異性散射效應(yīng)通常由一個簡單的輸運(yùn)修正截面來進(jìn)行各向同性近似,以提高計算效率。在多數(shù)常規(guī)的輕水堆物理計算中,如壓水堆的燃料組件計算,由于其泄漏較低,輸運(yùn)修正是充分合理的。但是,有多項研究表明[1-2],對于小堆或帶有控制棒的燃料組件,基于輸運(yùn)修正的中子輸運(yùn)計算的精度較差。典型問題包括Babcock &Wilcox (B&W)臨界實驗基準(zhǔn)題中的B&W1484[3],應(yīng)用傳統(tǒng)的輸運(yùn)修正造成的有效增殖因數(shù)偏差高達(dá)數(shù)千pcm。臨界實驗被廣泛應(yīng)用于中子輸運(yùn)程序的驗證與確認(rèn)工作中,是檢驗程序準(zhǔn)確性與有效性的重要手段[4-5],同時,臨界實驗裝置在中子學(xué)上通常具有強(qiáng)泄漏、強(qiáng)非均勻性的特征,在臨界實驗裝置的校核計算中必須考慮中子各向異性散射,而并非直接引用傳統(tǒng)的各向同性散射近似[6]。一些國際知名的燃料組件程序[6-8]在開發(fā)中引入了高階各向異性散射,然而引入高階散射,尤其是對于特征線方法(MOC)來說,意味著內(nèi)存和計算時間的大量增加。

針對這一問題,本文對高階各向異性散射的特征線輸運(yùn)計算方法開展研究,并實現(xiàn)異構(gòu)并行,選用LCT-011臨界實驗基準(zhǔn)進(jìn)行驗證計算。通過與蒙特卡羅計算程序結(jié)果的比較,確認(rèn)程序?qū)嶙V臨界實驗裝置的適用性與計算精度。此外,還進(jìn)行程序在高階散射計算中的性能分析,包括顯存占用以及計算效率。

1 計算方法

1.1 數(shù)值方法

MOC示意圖如圖1所示,在區(qū)域i內(nèi),沿著m方向的第k條特征線的中子輸運(yùn)方程為:

圖1 MOC示意圖Fig.1 Diagram of MOC

(1)

式中:φi,m,k(s)為沿著m方向的第k條特征線在s處的中子通量密度;Qi,m為總源項,包括裂變源與散射源;Σt,i為區(qū)域i的總宏觀截面。

(2)

進(jìn)而可推導(dǎo)出區(qū)域i的平均角通量密度[9]為:

(3)

式中,δAk為特征線段si,k的線寬。

各向異性源項的計算需要高階角通量密度矩的參與[9-10]:

(4)

其中:

(5)

μp=cosθp

(6)

(7)

對式(3)中的源項采用勒讓德展開:

(8)

本文僅考慮P1散射,則散射源項[9]可表示為:

(9)

(10)

1.2 程序開發(fā)

基于以上方法在已有的基于異構(gòu)平臺的特征線輸運(yùn)計算程序ALPHA上進(jìn)行開發(fā), ALPHA的開發(fā)情況參考文獻(xiàn)[11-14]。相比于P0散射計算,P1計算需要存儲各向異性散射源項以及平源區(qū)的平均角通量密度矩,因此平源區(qū)標(biāo)通量密度的計算方式從原先的在線原子操作變?yōu)閽呙柰戤吅蠼y(tǒng)一歸并。各向異性散射源項的計算與各項同性散射源計算方法差別不大,僅增加了角度循環(huán)與勒讓德展開階數(shù)循環(huán),整體計算流程示于圖2。

圖2 P1 MOC計算流程圖Fig.2 Computational flowchart of P1 MOC

2 臨界實驗描述

本文選用ICSBEP中的LEU-COMP-THERM-011(LCT011)臨界實驗基準(zhǔn)[15]進(jìn)行程序的驗證計算。

LCT011是1977—1979年在Babcock-Wilcox(B&W)林奇堡研究中心的CX-10臨界裝置上進(jìn)行的一系列使用鋁包殼以及低富集度氧化鈾的臨界實驗。LCT011中共有15組實驗,各組實驗使用相同富集度的燃料以及相同的結(jié)構(gòu)材料,但在活性區(qū)整體幾何形式、燃料棒陣列之間的柵距、B4C吸收棒位置以及慢化劑中硼酸濃度上體現(xiàn)差別。其中實驗1與實驗2采用相同的材料參數(shù),而堆芯幾何布置上分別采用Core Ⅰ和Core Ⅱ,如圖3所示。

圖3 LCT011臨界實驗Core Ⅰ(a)和Core Ⅱ(b)堆芯布置示意圖Fig.3 Loading diagram of Core Ⅰ (a) and Core Ⅱ(b) in LCT011 critical experiment

實驗3～9選用相同的堆芯布置,均為Core Ⅲ,而它們的差別在于慢化劑中不同的硼酸濃度。相較于實驗2所使用的堆芯布置Core Ⅱ,Core Ⅲ在原有3×3的燃料棒陣列之間加上了一個燃料柵元寬度的水隙。實驗10選用的堆芯布置Core Ⅳ與Core Ⅲ大體相同,但是在水隙中插入了一定數(shù)量的B4C吸收棒。Core Ⅲ與Core Ⅳ如圖4所示。

圖4 LCT011臨界實驗Core Ⅲ(a)和Core Ⅳ(b)堆芯布置示意圖Fig.4 Loading diagram of Core Ⅲ (a) and Core Ⅳ (b) in LCT011 critical experiment

實驗10～15選用的堆芯布置中燃料棒陣列之間的水隙逐漸增大,由1個燃料柵元寬度逐漸增大至4個燃料柵元寬度。在實驗11～14中相同堆芯布置的狀況下會改變B4C吸收棒的排布位置。各堆芯布置圖參考文獻(xiàn)[15]。燃料棒和B4C棒的軸向幾何結(jié)構(gòu)如圖5所示。

圖5 燃料棒(a)和B4C棒(b)幾何結(jié)構(gòu)示意圖Fig.5 Geometry diagram of fuel rod (a) and B4C rod (b)

綜上所述,LCT011中各臨界實驗具有很強(qiáng)的相關(guān)性,適合用于中子輸運(yùn)程序的驗證。

3 數(shù)值驗證

堆芯活性區(qū)徑向采用30 cm厚的反射層可以充分考慮對特征值和功率分布計算的影響,燃料柵元邊長為1.636 cm,因此,本文建模時在活性區(qū)外增加19個水柵元作為反射層,并將靠近活性區(qū)的10層?xùn)旁涂拷吔绲?層?xùn)旁x予不同的網(wǎng)格劃分。各柵元網(wǎng)格劃分方案列于表1。所有算例均在GPU服務(wù)器上完成,硬件配置列于表2。受限于硬件條件,本文選用2D 1/4堆芯進(jìn)行計算,最高散射階數(shù)取P1,并與OpenMC[16]計算結(jié)果進(jìn)行比較,程序與OpenMC的計算參數(shù)列于表3。

表1 LCT011計算中空間網(wǎng)格劃分方案Table 1 Spatial meshing scheme in LCT011 simulation

表2 計算平臺硬件配置及環(huán)境Table 2 Computing platform hardware configuration and environment

表3 OpenMC的計算參數(shù)Table 3 Computing parameter of OpenMC

3.1 keff計算結(jié)果

表4列出LCT011中15個算例的keff計算結(jié)果以及相對OpenMC參考解的誤差。相較于P0計算,P1計算可明顯提升計算精度,算例10的keff誤差最小,為-61 pcm;算例15的keff誤差最大,為378 pcm;15個算例平均keff誤差為236.5 pcm。整體來看,程序計算所得keff與OpenMC結(jié)果符合較好。

表4 LCT011的keff計算結(jié)果Table 4 keff calculation result of LCT011

3.2 棒功率計算結(jié)果

棒功率分布結(jié)果列于表5,未插入B4C棒的算例1～9以及算例15的平均棒功率相對誤差大都可以控制在0.5%以內(nèi),個別算例的平均相對誤差接近1%,除算例3的最大棒功率相對誤差為4.36%外,其他算例最大相對誤差均可控制在2%以內(nèi)。B4C棒的插入會增大計算誤差,算例10～13中最大棒功率相對誤差均大于10%,最大可至40.57%,而棒功率平均相對誤差均可控制在5%以內(nèi)。

表5 LCT011的棒功率誤差Table 5 Error of pin power of LCT011

圖6示出算例1、3、10、15的棒功率誤差分布,算例1和15代表著普通堆芯布置與大水隙布置,算例3與10分別為加入硼酸與插入B4C棒算例中棒功率誤差最大的算例。由圖6可知,與反射層或水隙相鄰位置棒功率誤差較大,而燃料陣列內(nèi)部棒功率相對誤差均在1%以內(nèi)。算例10中燃料陣列之間的水隙中插有B4C棒,而B4C棒所在位置燃料棒功率誤差較大,其他位置棒功率相對誤差均在5%以內(nèi)。在活性區(qū)反射層交界處以及B4C棒周圍因強(qiáng)泄漏和強(qiáng)吸收的原因,中子通量密度梯度大,導(dǎo)致這些位置棒功率誤差偏大。而這些位置的棒功率較低,且從圖6可看出,高通量密度梯度影響范圍具有局部性,因此整體功率以及keff的計算精度能維持在較高水平。

a——算例1;b——算例3;c——算例10;d——算例15圖6 LCT011的棒功率相對誤差分布Fig.6 Relative error distribution of pin power of LCT011

3.3 47群能譜計算結(jié)果

選用HELIOS 47群能群結(jié)構(gòu)[17],圖7示出算例1、3、10、15的能譜比較,該4個算例的全堆能譜平均相對誤差均在0.5%以內(nèi)。由圖7可知,共振群及快群能譜誤差較低且誤差隨能量分布平緩,均可控制在2%以內(nèi),且大都在0%附近。熱群及超熱群能譜誤差略大且有波動,但是這部分通量密度偏低,所以對整體計算結(jié)果影響較小。相較于無B4C棒的算例,含有B4C棒的算例在熱群與超熱群能譜誤差偏高,但在快群區(qū)間誤差與其他算例保持一致。在共振計算過程中,B4C引入了較大的各向異性,影響平源區(qū)宏觀截面歸并精度,因此在相應(yīng)能群處能譜誤差略微偏高,共振計算細(xì)節(jié)參考文獻(xiàn)[18]。

a——算例1;b——算例3;c——算例10;d——算例15圖7 LCT011的能譜比較結(jié)果Fig.7 Spectrum comparison result of LCT011

3.4 性能分析

在進(jìn)行P1計算時,源項計算與輸運(yùn)計算均在GPU上完成。為使得該性能分析可為三維計算服務(wù),計算時采用6極角進(jìn)行計算。各算例顯存使用情況與輸運(yùn)計算時間列于表6,所有算例的計算時間均可控制在7 min之內(nèi),具備較高的計算效率。表6列出總顯存消耗、各向異性散射源顯存消耗以及高階通量密度矩的顯存消耗。由表6可看出散射源項所占顯存占比很高,超過40%,散射源項的數(shù)據(jù)結(jié)構(gòu)與平源區(qū)角通量數(shù)據(jù)結(jié)構(gòu)相同,所以二者占用了總消耗顯存的80%以上。程序使用MOC-EX作為三維輸運(yùn)算法[19],該方法依靠軸向上具有幾何映射關(guān)系的平源區(qū)角通量密度相互傳遞進(jìn)行準(zhǔn)三維計算,所以平源區(qū)角通量密度所消耗的顯存不在未來優(yōu)化范圍目標(biāo)之內(nèi)。由此可得出結(jié)論,制約程序進(jìn)行更大規(guī)模臨界實驗裝置物理計算的瓶頸在于散射源項的存儲。各向異性散射源項可以放置在輸運(yùn)掃描過程中在線計算,但該方法對計算效率的損害需要進(jìn)一步分析。

4 結(jié)論

針對高階各向異性散射的特征線方法開展了研究,并在異構(gòu)計算平臺上進(jìn)行了程序?qū)崿F(xiàn)。基于LCT011臨界實驗對程序的高階散射輸運(yùn)計算功能進(jìn)行了驗證及分析。數(shù)值結(jié)果表明,該方法和程序計算熱譜臨界實驗裝置有效增殖因數(shù)與蒙特卡羅所得結(jié)果較為吻合。在無B4C棒插入的情況下計算所得棒功率平均相對誤差可控制在0.5%以內(nèi),有B4C棒的情況下平均相對誤差可控制在5%以內(nèi)。從能譜結(jié)果來看,無論B4C插入與否,程序計算所得能譜平均相對誤差可控制在0.5%以內(nèi)。計算LCT011的性能分析表明,程序能夠在較短時間內(nèi)完成二維臨界裝置的精細(xì)物理計算,但各向異性散射源項是制約程序進(jìn)一步擴(kuò)大計算規(guī)模的瓶頸。綜上所述,程序在計算熱譜臨界裝置時具有與蒙特卡羅程序同等精度,且具有較高計算效率。未來工作將聚焦于高階散射計算條件下的異構(gòu)計算平臺上的顯存優(yōu)化以及效率優(yōu)化。