張顯 劉仕倡 魏軍俠 李樹 王鑫3) 上官丹驊

1) (北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100094)

2) (華北電力大學(xué)核科學(xué)與工程學(xué)院,北京 102206)

3) (中國(guó)工程物理研究院高性能數(shù)值模擬軟件中心,北京 100088)

全局計(jì)數(shù)問(wèn)題在反應(yīng)堆pin-by-pin 模型蒙特卡羅模擬和多物理耦合計(jì)算中動(dòng)態(tài)粒子輸運(yùn)蒙特卡羅模擬等重大研究領(lǐng)域中都有廣泛的應(yīng)用場(chǎng)景.大量的全局減方差算法研究立足于全局計(jì)數(shù)誤差分布的展平,由此提高全局計(jì)數(shù)的整體效率.本工作針對(duì)兩種高效全局減方差算法,即均勻計(jì)數(shù)密度算法(屬于源偏倚算法的一種)和權(quán)窗算法的結(jié)合展開研究,提出利用均勻計(jì)數(shù)密度算法的偏倚因子調(diào)整權(quán)窗下限,由此實(shí)現(xiàn)兩種算法的有機(jī)結(jié)合.基于Hoogenboom-Martin 壓水堆全堆基準(zhǔn)題中開展了一系列對(duì)比測(cè)試,驗(yàn)證了混合全局減方差算法更優(yōu)于單一權(quán)窗算法或均勻計(jì)數(shù)密度算法,尤其是在降低最大誤差方面.同時(shí),基于新的指標(biāo),驗(yàn)證了均勻計(jì)數(shù)密度算法較經(jīng)典的均勻裂變?cè)此惴ň哂懈玫谋憩F(xiàn).研究結(jié)果表明,本文提出的混合全局減方差算法能高效求解全局計(jì)數(shù)問(wèn)題,進(jìn)一步促進(jìn)了相關(guān)領(lǐng)域的研究.

1 引言

蒙特卡羅(Monte Carlo,MC)方法具有幾何建模能力強(qiáng)、物理過(guò)程描述高保真等優(yōu)點(diǎn),被廣泛應(yīng)用于定態(tài)和動(dòng)態(tài)粒子輸運(yùn)問(wèn)題的模擬.隨著研究的深入,輸運(yùn)問(wèn)題的幾何模型愈加精細(xì),考慮的因素越來(lái)越多,例如反應(yīng)堆pin-by-pin 模型[1,2]和多物理耦合計(jì)算中大規(guī)模動(dòng)態(tài)粒子輸運(yùn)模型等[3].龐大的幾何和計(jì)數(shù)規(guī)模以及高效高精度的計(jì)算需求,給MC 粒子輸運(yùn)模擬帶來(lái)巨大挑戰(zhàn).由于所模擬系統(tǒng)的空間不均勻性,計(jì)數(shù)統(tǒng)計(jì)誤差在全局范圍內(nèi)呈現(xiàn)不均勻分布,由此帶來(lái)整體效率的低下.解決這一問(wèn)題的本質(zhì)困難在于高功率區(qū)域因具有較多的粒子樣本數(shù),能較快獲得統(tǒng)計(jì)收斂的結(jié)果,而低功率區(qū)域的收斂則耗時(shí)巨大,僅單純?cè)黾訕颖究倲?shù),只會(huì)導(dǎo)致絕大部分計(jì)算資源浪費(fèi)在已收斂的高功率區(qū)域,難以(在有些情況下是不可能的)獲得全局收斂的計(jì)數(shù)結(jié)果,而這些結(jié)果對(duì)于反應(yīng)堆計(jì)算或高置信度多物理耦合計(jì)算至關(guān)重要.

為實(shí)現(xiàn)全局計(jì)數(shù)整體誤差分布的展平,需要引入相關(guān)的全局減方差算法指導(dǎo)MC 粒子輸運(yùn),提高全局計(jì)數(shù)的整體效率[4,5].均勻裂變?cè)?UFS)算法[6–8]是針對(duì)臨界計(jì)算提出的一種高效全局減方差算法,根據(jù)裂變中子源的密度分布重新分配裂變?cè)?以便在低功率區(qū)域人為地產(chǎn)生更多的裂變中子.基于UFS 算法的啟發(fā),均勻計(jì)數(shù)密度(UTD)算法[9]被提出,其利用目標(biāo)計(jì)數(shù)密度來(lái)指導(dǎo)源粒子的偏倚,獲得了更高的全局減方差性能.此外,權(quán)窗(WW)算法[10,11]也是一種被廣泛應(yīng)用的全局減方差算法,不同于UTD 和UFS 算法,WW 算法是對(duì)粒子的輸運(yùn)過(guò)程進(jìn)行偏倚,以引導(dǎo)粒子輸運(yùn)到更廣泛的區(qū)域.

為進(jìn)一步提高M(jìn)C 臨界計(jì)算全局計(jì)數(shù)問(wèn)題的整體效率,本文提出一種結(jié)合UTD 算法和WW算法的混合全局減方差方法,其利用UTD 的偏倚因子動(dòng)態(tài)調(diào)整WW 下限,利用WW 減小UTD 方法引起的權(quán)重波動(dòng),以此實(shí)現(xiàn)兩種算法的有機(jī)結(jié)合.這一方法在MC 粒子輸運(yùn)程序cosRMC[12–14]上進(jìn)行了驗(yàn)證.第2 節(jié)介紹了UFS 和UTD 算法的基本思想;第3 節(jié)對(duì)混合算法的思想和實(shí)現(xiàn)方法進(jìn)行了描述;第4 節(jié)基于新的指標(biāo)深入研究了UTD算法的效率,并開展混合算法的測(cè)試驗(yàn)證;第5 節(jié)給出了相應(yīng)結(jié)論.

2 UFS 和UTD 算法原理

在反應(yīng)堆模擬計(jì)算中,不同幾何柵元間的功率密度會(huì)有較大差異,全堆中子樣本數(shù)量就會(huì)呈現(xiàn)不均勻分布,導(dǎo)致全局計(jì)數(shù)不能同步收斂.UFS 算法的基本思想: 在保證結(jié)果無(wú)偏的前提下,對(duì)裂變中子源分布進(jìn)行調(diào)整.由于MC 方法進(jìn)行臨界計(jì)算是以迭代形式開展的,上一代產(chǎn)生的次級(jí)裂變中子將作為下一代的裂變中子源.基于這一特點(diǎn),UFS算法根據(jù)當(dāng)前代的裂變中子數(shù)密度分布產(chǎn)生偏倚因子,在下一代開始時(shí)指導(dǎo)裂變?cè)捶植嫉恼{(diào)整.為便于描述裂變中子源分布,在堆芯區(qū)域疊加均勻網(wǎng)格對(duì)空間進(jìn)行離散,以網(wǎng)格為單元執(zhí)行源粒子的偏倚.UFS 偏倚因子的設(shè)置方法為

其中,Nt為總裂變?cè)粗凶訑?shù),m為總網(wǎng)格數(shù),Ni為網(wǎng)格i內(nèi)的裂變?cè)粗凶訑?shù).

引入源偏倚因子βi后,網(wǎng)格i內(nèi)每次碰撞產(chǎn)生的裂變中子數(shù)[15]將被調(diào)整為

其中,w是發(fā)生碰撞的中子權(quán)重;υ為平均次級(jí)裂變中子數(shù);為宏觀裂變中子截面;是宏觀總截面.為保證計(jì)算結(jié)果無(wú)偏,下一代裂變?cè)粗凶拥臋?quán)重ws將調(diào)整為ws/βi.

上述算法將導(dǎo)致低功率區(qū)域分裂出更多的小權(quán)重中子,而高功率區(qū)域則相應(yīng)減少了裂變中子數(shù),因此不會(huì)增加額外的計(jì)算耗時(shí).如果減方差目標(biāo)是展平某種全局計(jì)數(shù)的統(tǒng)計(jì)誤差分布,以目標(biāo)計(jì)數(shù)密度指導(dǎo)源粒子的偏倚可能比基于裂變中子數(shù)密度的偏倚效率更高[9].基于此,UTD 算法提出偏倚因子的設(shè)定方法為

其中Tt為所有網(wǎng)格目標(biāo)計(jì)數(shù)值之和;m為總網(wǎng)格數(shù);Ti為網(wǎng)格i的目標(biāo)計(jì)數(shù)值.上述兩種算法本質(zhì)上都是源偏倚算法.

3 基于UTD 和WW 的混合全局減方差算法

WW 算法是一種基于分裂和輪盤賭的全局減方差方法,也需要借助網(wǎng)格來(lái)為不同空間區(qū)域提供WW.每一個(gè)網(wǎng)格的WW 由三個(gè)參數(shù)組成,包括WW 上限、WW 下限和輪盤賭存活權(quán)重.每當(dāng)粒子到達(dá)柵元邊界、碰撞點(diǎn)以及飛行每個(gè)平均自由程后,都會(huì)對(duì)粒子的權(quán)重進(jìn)行檢查.如圖1 所示,如果粒子權(quán)重低于WW 下限,就會(huì)觸發(fā)輪盤賭機(jī)制,有效地截?cái)嘈?quán)重的粒子;如果粒子權(quán)重高于WW 上限,對(duì)粒子執(zhí)行分裂操作,增加粒子樣本數(shù).通過(guò)為低功率區(qū)域設(shè)置較小的WW 參數(shù),為高功率區(qū)域設(shè)置較大的WW 參數(shù),可以實(shí)現(xiàn)計(jì)算資源的均勻分配.WW 算法是一種輸運(yùn)過(guò)程偏倚算法.

圖1 權(quán)窗原理Fig.1.Working principle of weight window.

為結(jié)合源偏倚與輸運(yùn)過(guò)程偏倚各自的優(yōu)勢(shì),獲得臨界計(jì)算全局計(jì)數(shù)整體效率的進(jìn)一步提高,本文提出一種基于UTD 和WW 的混合算法.由于UTD 算法會(huì)改變裂變中子的權(quán)重,可能會(huì)引起較大的粒子權(quán)重波動(dòng),不利于統(tǒng)計(jì)結(jié)果的整體收斂,而WW 算法可以將粒子權(quán)重控制在合理范圍內(nèi),因此混合算法預(yù)計(jì)可以獲得更高的整體效率.

UTD 方法和WW 方法均基于網(wǎng)格執(zhí)行偏倚操作,因此在混合算法中兩者可以共用一套網(wǎng)格劃分方案.在低功率區(qū)域,UTD 和WW 方法都會(huì)分裂中子,混合算法將建立兩個(gè)臨時(shí)儲(chǔ)存庫(kù),對(duì)這些粒子進(jìn)行臨時(shí)分類存放,按序完成所有粒子的輸運(yùn)模擬.由于UTD 算法在低功率區(qū)域會(huì)分裂出極小權(quán)重的中子,WW 的輪盤賭機(jī)制可能直接截?cái)噙@類粒子,對(duì)UTD 算法的效果造成一定削弱.因此,提出使用UTD 偏倚因子βi來(lái)調(diào)整網(wǎng)格WW 下限WL:

通過(guò)這種方法,在不同功率區(qū)域根據(jù)UTD 偏倚因子,合理地降低或抬高WW 下限,可以減少WW對(duì)UTD 性能的負(fù)面影響,實(shí)現(xiàn)兩種方法的有機(jī)結(jié)合.

4 數(shù)值結(jié)果及分析

選擇在Hoogenboom-Martin 壓水堆全堆基準(zhǔn)題[16,17]上開展相關(guān)的測(cè)試驗(yàn)證.如圖2 所示,該模型堆芯徑向半徑為209 cm,軸向高度為366 cm,共包含241 個(gè)燃料組件,燃料組件為 17×17 布置;每個(gè)組件內(nèi)呈現(xiàn)17×17 的棒分布,包含264 個(gè)燃料棒和25 個(gè)控制棒通道.從圖2(a)和圖2(b)(不采用任何全局減方差算法)可以看出,基準(zhǔn)模型的徑向功率分布具有顯著不均勻性,導(dǎo)致統(tǒng)計(jì)誤差分布也呈現(xiàn)嚴(yán)重不均.將堆芯沿橫向和縱向劃分成289×289 的網(wǎng)格,其中燃料區(qū)網(wǎng)格共計(jì)69649 個(gè).計(jì)算條件為非活躍代數(shù)200,活躍代數(shù)300,每代初始粒子數(shù)50000,采用50 核并行計(jì)算,統(tǒng)計(jì)每個(gè)網(wǎng)格的中子裂變功率.

圖2 Hoogenboom-Martin 基準(zhǔn)題(a) 幾何橫截面;(b) 功率分布(MW);(c) 統(tǒng)計(jì)誤差分布Fig.2.Hoogenboom Martin benchmark: (a) Geometric cross-section;(b) power distribution (MW);(c) statistical error distribution.

4.1 基于新指標(biāo)的UTD 和UFS 對(duì)比分析

為更清晰地了解UTD 算法的優(yōu)勢(shì)所在,基于新指標(biāo)對(duì)UTD 和UFS 進(jìn)行了對(duì)比分析.在H-M基準(zhǔn)題計(jì)算中,堆芯模型具有1/4 對(duì)稱性,4 個(gè)對(duì)稱區(qū)域的物理量在理想情況下應(yīng)完全相同,而由于統(tǒng)計(jì)誤差的存在,導(dǎo)致堆芯物理對(duì)稱區(qū)域的計(jì)算結(jié)果略有不同,稱之為計(jì)算不對(duì)稱性.引入變異系數(shù)Cυ定量描述這種計(jì)算不對(duì)稱程度,變異系數(shù)越大表示不對(duì)稱性程度越大[18]:

其中,S和分別是對(duì)稱區(qū)域4 個(gè)計(jì)數(shù)量的標(biāo)準(zhǔn)偏差和平均.

圖3 給出了Basic,UFS 和UTD 三種情況下的變異系數(shù)分布.從圖3 可以看出,不使用任何源偏倚的Basic 情況下,堆芯外圍的計(jì)算不對(duì)稱程度明顯比中心的不對(duì)稱程度大;使用UFS 算法時(shí),堆芯外圍的Cυ明顯降低,計(jì)算不對(duì)稱程度較Basic的小;UTD 算法下的計(jì)算不對(duì)稱度相較UFS 有了更進(jìn)一步的改善.

圖3 變異系數(shù)分布(a) Basic;(b) UFS;(c) UTDFig.3.Distribution of the coefficient of variation: (a) Basic;(b) UFS;(c) UTD.

此外,UFS 算法和UTD 算法在臨界計(jì)算的每個(gè)活躍代都會(huì)對(duì)偏倚因子進(jìn)行更新,選取了堆芯中橫向289 個(gè)連續(xù)的網(wǎng)格,計(jì)算得到UTD 和UFS算法在這些網(wǎng)格中的偏倚因子的方差分布,見(jiàn)圖4.UTD 偏倚因子的方差整體小于UFS 偏倚因子的方差,UTD 算法的偏倚因子在迭代過(guò)程中波動(dòng)更小,表明UTD 算法相比UFS 算法更具穩(wěn)定性.

圖4 偏倚因子的方差分布Fig.4.Variance distribution of bias factors.

為了進(jìn)一步對(duì)比兩者的減方差效果,引入品質(zhì)因子FOM_MAX 和FOM_95 來(lái)表征計(jì)算效率[19,20]:

這里,T是計(jì)算時(shí)間,Remax是所有網(wǎng)格計(jì)數(shù)的統(tǒng)計(jì)誤差的最大值,Re95表示某一網(wǎng)格計(jì)數(shù)的統(tǒng)計(jì)誤差,其使至少95%的網(wǎng)格計(jì)數(shù)的統(tǒng)計(jì)誤差都不大于該值.表1 給出了三種計(jì)算條件下中子注量率的統(tǒng)計(jì)誤差和品質(zhì)因子.從兩種品質(zhì)因子來(lái)看,UTD 算法和UFS 算法均獲得了計(jì)算效率的提高,并且UTD 算法的計(jì)算效率略優(yōu)于UFS 算法.在降低最大誤差的問(wèn)題上,UTD 算法的計(jì)算效率是UFS 算法的1.36 倍.

表1 UTD 算法和UFS 算法的計(jì)算結(jié)果對(duì)比Table 1.Comparison of calculation results of UTD and UFS.

4.2 混合算法的驗(yàn)證

統(tǒng)計(jì)誤差的累積分布如圖5 所示,UTD 算法下的統(tǒng)計(jì)誤差較集中落在3.7%—5.0%區(qū)間內(nèi),混合算法和WW 算法下的統(tǒng)計(jì)誤差較集中落在2.7%—3.2%區(qū)間內(nèi).從結(jié)果可以看出,WW 算法下的統(tǒng)計(jì)誤差整體小于UTD 算法,且最大和最小誤差的差值相比UTD 算法也更小,說(shuō)明WW 算法的減方差力度比UTD 算法更大.使用混合算法時(shí),UTD 算法會(huì)在臨界計(jì)算中的每個(gè)活躍代開始對(duì)裂變?cè)捶植歼M(jìn)行調(diào)整,WW 算法會(huì)在粒子輸運(yùn)過(guò)程中對(duì)粒子進(jìn)行偏倚,兩者的共同作用使得低功率區(qū)域具有更多的粒子樣本,在低功率區(qū)域?qū)崿F(xiàn)更進(jìn)一步的減方差效果,從而混合算法的最大統(tǒng)計(jì)誤差小于WW 算法和UTD 算法.表2 列出了三種減方差算法的整體效率,通過(guò)品質(zhì)因子可以得出,針對(duì)降低最大誤差的問(wèn)題,混合算法的計(jì)算效率分別是WW 算法和UTD 算法的2.6 倍和3 倍.

表2 計(jì)算結(jié)果對(duì)比Table 2.Comparison of calculation results.

圖5 統(tǒng)計(jì)誤差的累積分布Fig.5.Cumulative distribution of statistical errors.

5 結(jié)論

圍繞MC 粒子輸運(yùn)模擬全局計(jì)數(shù)問(wèn)題統(tǒng)計(jì)誤差分布不均的問(wèn)題,本文提出了一種結(jié)合均勻計(jì)數(shù)密度算法和WW 算法的混合算法,通過(guò)引入WW減少了均勻計(jì)數(shù)密度算法導(dǎo)致的權(quán)重波動(dòng),其WW 下限利用均勻計(jì)數(shù)密度算法偏倚因子進(jìn)行調(diào)整,本質(zhì)上實(shí)現(xiàn)了源偏倚和輸運(yùn)過(guò)程偏倚的有機(jī)結(jié)合.在Hoogenboom-Martin 基準(zhǔn)題的驗(yàn)證計(jì)算中,基于新的指標(biāo)對(duì)比分析了均勻計(jì)數(shù)密度算法和UFS 算法,進(jìn)一步驗(yàn)證了UTD 算法的高效性.同時(shí),計(jì)算結(jié)果表明,混合算法的整體效率較均勻計(jì)數(shù)密度算法或WW 算法有進(jìn)一步的提高.在降低最大誤差方面,混合算法的整體效率分別是WW算法和均勻計(jì)數(shù)密度算法的2.6 倍和3倍,驗(yàn)證了混合算法的優(yōu)越性.