異構并行的高階散射特征線方法及其在臨界實驗裝置模擬中的應用

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

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

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

針對這一問題,本文對高階各向異性散射的特征線輸運計算方法開展研究,并實現異構并行,選用LCT-011臨界實驗基準進行驗證計算。通過與蒙特卡羅計算程序結果的比較,確認程序對熱譜臨界實驗裝置的適用性與計算精度。此外,還進行程序在高階散射計算中的性能分析,包括顯存占用以及計算效率。

1 計算方法

1.1 數值方法

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

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

(1)

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

(2)

進而可推導出區域i的平均角通量密度[9]為:

(3)

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

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

(4)

其中:

(5)

μp=cosθp

(6)

(7)

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

(8)

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

(9)

(10)

1.2 程序開發

基于以上方法在已有的基于異構平臺的特征線輸運計算程序ALPHA上進行開發, ALPHA的開發情況參考文獻[11-14]。相比于P0散射計算,P1計算需要存儲各向異性散射源項以及平源區的平均角通量密度矩,因此平源區標通量密度的計算方式從原先的在線原子操作變為掃描完畢后統一歸并。各向異性散射源項的計算與各項同性散射源計算方法差別不大,僅增加了角度循環與勒讓德展開階數循環,整體計算流程示于圖2。

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

2 臨界實驗描述

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

LCT011是1977—1979年在Babcock-Wilcox(B&W)林奇堡研究中心的CX-10臨界裝置上進行的一系列使用鋁包殼以及低富集度氧化鈾的臨界實驗。LCT011中共有15組實驗,各組實驗使用相同富集度的燃料以及相同的結構材料,但在活性區整體幾何形式、燃料棒陣列之間的柵距、B4C吸收棒位置以及慢化劑中硼酸濃度上體現差別。其中實驗1與實驗2采用相同的材料參數,而堆芯幾何布置上分別采用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 Ⅲ大體相同,但是在水隙中插入了一定數量的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吸收棒的排布位置。各堆芯布置圖參考文獻[15]。燃料棒和B4C棒的軸向幾何結構如圖5所示。

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

綜上所述,LCT011中各臨界實驗具有很強的相關性,適合用于中子輸運程序的驗證。

3 數值驗證

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

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

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

表3 OpenMC的計算參數Table 3 Computing parameter of OpenMC

3.1 keff計算結果

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

表4 LCT011的keff計算結果Table 4 keff calculation result of LCT011

3.2 棒功率計算結果

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

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

圖6示出算例1、3、10、15的棒功率誤差分布,算例1和15代表著普通堆芯布置與大水隙布置,算例3與10分別為加入硼酸與插入B4C棒算例中棒功率誤差最大的算例。由圖6可知,與反射層或水隙相鄰位置棒功率誤差較大,而燃料陣列內部棒功率相對誤差均在1%以內。算例10中燃料陣列之間的水隙中插有B4C棒,而B4C棒所在位置燃料棒功率誤差較大,其他位置棒功率相對誤差均在5%以內。在活性區反射層交界處以及B4C棒周圍因強泄漏和強吸收的原因,中子通量密度梯度大,導致這些位置棒功率誤差偏大。而這些位置的棒功率較低,且從圖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群能譜計算結果

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

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

3.4 性能分析

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

4 結論

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