王學(xué)棟,郭丁情,2,佟立麗,曹學(xué)武,*

(1.上海交通大學(xué) 機(jī)械與動力工程學(xué)院,上海200240;2.中廣核工程有限公司,廣東 深圳518000)

核電廠嚴(yán)重事故下輻射環(huán)境是評估電廠儀表、設(shè)備可用性和人員可達(dá)性的重要參考,目前的嚴(yán)重事故管理導(dǎo)則及核電廠安全規(guī)程都要求,在可能的范圍內(nèi),重要儀表和設(shè)備能夠以合理的可信度在嚴(yán)重事故條件下達(dá)到設(shè)計(jì)要求,以及評估緩解事故的操作人員可能到達(dá)的電廠區(qū)域,以確定其環(huán)境條件(高溫、高壓、高輻射劑量等)能夠支持緩解措施的執(zhí)行[1,2]。因此,分析核電廠在嚴(yán)重事故下的輻射環(huán)境、評估重要區(qū)域的劑量率水平對于核電廠嚴(yán)重事故管理、應(yīng)急響應(yīng)等環(huán)節(jié)具有重要意義。當(dāng)前核電廠屏蔽結(jié)構(gòu)的復(fù)雜性要求輻射屏蔽計(jì)算方法能夠在較為逼真地還原實(shí)際模型的情況下,盡量減少因結(jié)構(gòu)簡化帶來的計(jì)算偏差,且模擬結(jié)果能夠以合理的置信水平指導(dǎo)相關(guān)的技術(shù)問題,此外,方法的選取與應(yīng)用也必須考慮到計(jì)算時(shí)間成本的維護(hù),盡可能地避免計(jì)算過程的冗長耗時(shí)給輻射屏蔽問題本身帶來的負(fù)面影響。基于上述客觀要求,蒙特卡羅方法(MC)是當(dāng)前屏蔽設(shè)計(jì)與分析領(lǐng)域解決復(fù)雜幾何結(jié)構(gòu)建模問題的首選工具。

MC方法的有效應(yīng)用要求其在保持固有優(yōu)勢的基礎(chǔ)上,盡可能地提高收斂速度。復(fù)雜核能系統(tǒng)的屏蔽計(jì)算首先面臨的是粒子輸運(yùn)過程中的深穿透問題,若是直接采用MC方法,在有限的計(jì)算資源和時(shí)間成本控制下,這種深穿透特點(diǎn)會給目標(biāo)區(qū)域探測器的計(jì)數(shù)結(jié)果造成很大的統(tǒng)計(jì)誤差,往往使得模擬結(jié)果不可信,甚至出現(xiàn)計(jì)數(shù)為零的極端情況。在這種形勢下,采取合理的減方差方法以使得模擬結(jié)果的統(tǒng)計(jì)誤差處于可接受的范圍逐漸成為了當(dāng)前MC方法應(yīng)用研究的重點(diǎn)課題。傳統(tǒng)的減方差技巧的應(yīng)用將落腳點(diǎn)放在了屏蔽系統(tǒng)的局部區(qū)域,通常只對局部區(qū)域計(jì)算結(jié)果統(tǒng)計(jì)精度的提高有效,而當(dāng)前核能系統(tǒng)的屏蔽計(jì)算往往會涉及到全局范圍,此外,相關(guān)參數(shù)的設(shè)置通常會依賴于程序使用者的模擬計(jì)算經(jīng)驗(yàn),缺乏有效的理論支撐,導(dǎo)致這些技巧的應(yīng)用效果十分有限。因此,如何從全局角度提高M(jìn)C方法應(yīng)用的準(zhǔn)確性和時(shí)效性是亟需探索和解決的關(guān)鍵問題。普遍認(rèn)為基于粒子分裂和輪盤賭技巧的輸運(yùn)偏倚是提高蒙特卡羅方法全局計(jì)算精度和效率的最佳手段[3,4]。利用MCNP程序的權(quán)重窗功能可對粒子在輸運(yùn)過程中進(jìn)行一定程度的偏倚,通過控制粒子的權(quán)重,減少高計(jì)數(shù)區(qū)域的抽樣次數(shù),而增加低計(jì)數(shù)區(qū)域的抽樣次數(shù),從而實(shí)現(xiàn)展平蒙特卡羅模擬粒子在全局范圍內(nèi)的分布、降低計(jì)數(shù)結(jié)果統(tǒng)計(jì)誤差的目的[5]。應(yīng)用于權(quán)重窗計(jì)算的理論方法主要有基于確定性計(jì)算結(jié)果的伴隨函數(shù)方法[6-9]和基于蒙特卡羅正算結(jié)果的迭代方法[10]。伴隨理論的難點(diǎn)在于源偏倚和輸運(yùn)偏倚之間的耦合問題,以及前后兩次計(jì)算存在的網(wǎng)格匹配性問題,因而限制了其適用的問題類型范圍。蒙特卡羅正算迭代方法盡管存在誤差傳遞等不利因素,但因其實(shí)現(xiàn)手段較為簡單,且能有效提高蒙特卡羅程序最終計(jì)算結(jié)果的精度和品質(zhì)因數(shù)(FOM),目前已獲得廣泛的應(yīng)用。

本文以某重水堆核電廠在嚴(yán)重事故下的輻射環(huán)境為模擬對象,針對模擬結(jié)果收斂速度慢的缺點(diǎn),利用蒙特卡羅正算輸運(yùn)偏倚理論,通過多次迭代計(jì)算獲得基于網(wǎng)格通量計(jì)數(shù)的權(quán)重窗下限參數(shù),將其運(yùn)用到最終的蒙特卡羅模擬過程,探索出了一條在輻射屏蔽計(jì)算中運(yùn)用全局減方差方法以提高蒙特卡羅計(jì)算效果的途徑。

1 相關(guān)概念和理論

1.1 蒙特卡羅方法的統(tǒng)計(jì)不確定性

應(yīng)用蒙特卡羅方法得到的模擬結(jié)果本質(zhì)上是一個(gè)數(shù)學(xué)統(tǒng)計(jì)量,因此模擬結(jié)果必然存在著統(tǒng)計(jì)不確定性的問題。通常,蒙特卡羅模擬結(jié)果的統(tǒng)計(jì)不確定性用相對誤差估計(jì)R來表征,即68%置信水平下的估計(jì)值,用公式(1)和公式(2)式表示:

式中:Sˉx——樣本均值ˉx的標(biāo)準(zhǔn)差;

S2——樣本方差;

N——模擬粒子歷史數(shù)。

在中心極限定理適用的情況下,相對誤差估計(jì)的平方R2應(yīng)當(dāng)正比于1/N。由于蒙特卡羅計(jì)算過程的耗時(shí)T一般是與模擬的粒子歷史數(shù)成正相關(guān)的,因此,一個(gè)較好的蒙特卡羅模擬結(jié)果應(yīng)當(dāng)對應(yīng)著R2T趨于一個(gè)常數(shù)。為了定量化蒙特卡羅模擬過程的計(jì)算效率,文獻(xiàn)[11]提出了一個(gè)品質(zhì)因數(shù)(FOM)的概念,用公式(3)表示:

品質(zhì)因數(shù)的概念反映在計(jì)算效率上是其值與計(jì)算耗時(shí)之間的關(guān)系,可以作為一個(gè)表征參數(shù),用來量化和比較蒙特卡羅模擬計(jì)算中減方差方法的實(shí)際應(yīng)用效果。

1.2 全局減方差理論

對于深穿透問題,直接蒙特卡羅模擬無法獲得優(yōu)異的全局特性,一個(gè)重要表現(xiàn)就是全局模擬結(jié)果的統(tǒng)計(jì)不確定性分布很不均勻,深穿透區(qū)域由于抽樣數(shù)目很少甚至為零,導(dǎo)致其統(tǒng)計(jì)誤差很大,而在源項(xiàng)區(qū)域或近源區(qū),由于計(jì)數(shù)單元能夠獲得足夠的樣本數(shù)目,統(tǒng)計(jì)誤差往往很低甚至?xí)h(yuǎn)遠(yuǎn)超過計(jì)算要求,這種統(tǒng)計(jì)誤差分布的極度不均勻性在程序運(yùn)行層面來說是一種計(jì)算資源在問題系統(tǒng)不同區(qū)域分配極度不均的體現(xiàn),往往也會造成模擬計(jì)算耗時(shí)嚴(yán)重過長,大大增加了計(jì)算成本。因此,全局減方差的首要目標(biāo)就是使用一切可能的手段展平問題系統(tǒng)內(nèi)計(jì)算結(jié)果的統(tǒng)計(jì)誤差,提高深穿透區(qū)域的抽樣數(shù)目,相應(yīng)地減少源項(xiàng)區(qū)域和近源區(qū)的抽樣數(shù)目,具體在粒子輸運(yùn)過程的層面,是要采取合理的偏倚方法引導(dǎo)模擬粒子盡可能地向深穿透區(qū)域進(jìn)行輸運(yùn),即展平蒙特卡羅模擬粒子在全局范圍內(nèi)的分布,均勻分配計(jì)算資源。

假定全局問題的蒙特卡羅計(jì)算在使用網(wǎng)格權(quán)重窗后,網(wǎng)格柵元內(nèi)蒙特卡羅模擬粒子的平均權(quán)重與網(wǎng)格權(quán)重窗上限和下限參數(shù)的平均值成正相關(guān),由于在MCNP內(nèi)權(quán)重窗上限(Wsplit)一般設(shè)置為下限的若干倍,因此可構(gòu)建如公式(4)所示的權(quán)重窗下限參數(shù)的計(jì)算方法。

式中:Wth,i——網(wǎng)格i內(nèi)權(quán)重窗下限參數(shù);

φi——網(wǎng)格i內(nèi)的粒子通量計(jì)數(shù)結(jié)果;

λ——權(quán)重窗上限和下限參數(shù)之比,即Wsplit/Wth;

max(φ)——系統(tǒng)全局范圍內(nèi)各個(gè)計(jì)數(shù)網(wǎng)格最大的粒子通量。

由公式(4)可以看出,計(jì)數(shù)結(jié)果很小的網(wǎng)格柵元對應(yīng)的權(quán)重窗參數(shù)也很小,即粒子發(fā)生分裂的閾值很低,這樣就會促使更多的粒子發(fā)生分裂,增加?xùn)旁獌?nèi)的蒙特卡羅模擬粒子數(shù),提高這些網(wǎng)格內(nèi)的模擬計(jì)算時(shí)長,從而降低了計(jì)數(shù)結(jié)果的統(tǒng)計(jì)不確定性。相應(yīng)地,計(jì)數(shù)結(jié)果很大的網(wǎng)格柵元發(fā)生粒子分裂的閾值很高,限制了其模擬粒子數(shù)目的進(jìn)一步增加,同時(shí)由于發(fā)生輪盤賭的概率較大,會有相當(dāng)一部分低權(quán)重粒子的輸運(yùn)歷史停止,有效縮短了模擬運(yùn)行時(shí)長。綜合來看,利用網(wǎng)格通量來構(gòu)建權(quán)重窗參數(shù)能夠起到展平問題系統(tǒng)內(nèi)蒙特卡羅模擬粒子數(shù)目和均勻分配計(jì)算資源的作用。

2 模擬計(jì)算與方法應(yīng)用

2.1 計(jì)算模型

針對輻射環(huán)境的蒙特卡羅模擬計(jì)算,本文對某重水堆核電廠安全殼廠房、核輔助廠房和乏燃料廠房等所有對放射性光子屏蔽有重要貢獻(xiàn)的建筑設(shè)施進(jìn)行了空間幾何描述,幾何結(jié)構(gòu)涵蓋了廠房建筑的4個(gè)標(biāo)高層:93.9 m標(biāo)高層包括應(yīng)急堆芯冷卻系統(tǒng)(ECCS)泵、ECCS管道、應(yīng)急閘門和通道等建筑結(jié)構(gòu);100 m標(biāo)高層包括ECCS換熱房間和管道、走廊等結(jié)構(gòu)設(shè)施;105.41 m標(biāo)高層包括設(shè)備閘門等結(jié)構(gòu)設(shè)施;109.22 m標(biāo)高層包括主控室等房間。在建模過程中,忽略了不重要的附屬結(jié)構(gòu),并且在一定程度上對廠房結(jié)構(gòu)進(jìn)行了等效處理。模型的整體尺度為150 m×73 m×35 m,圖1給出了核電廠蒙特卡羅建模結(jié)果。

圖1 重水堆核電廠MC模型Fig.1 The MC model of heavy water reactornuclear power

源項(xiàng)模型選取為假想的嚴(yán)重事故下的輻射情景,假定堆芯釋放出的大量放射性核素都包容在安全殼內(nèi)沒有發(fā)生泄漏,且核素在安全殼內(nèi)均勻分布。源項(xiàng)釋放份額按照NUREG-1465給出的相關(guān)要求[12]進(jìn)行選取,以堆芯燃料棒開始發(fā)生裸露后0.5 h的時(shí)刻作為本文的計(jì)算基準(zhǔn)。

2.2 計(jì)算結(jié)果與討論

針對上述核電廠嚴(yán)重事故下輻射環(huán)境的計(jì)算模型,本文對其全廠范圍內(nèi)Gamma通量的空間分布進(jìn)行了蒙特卡羅模擬,劃分的計(jì)數(shù)網(wǎng)格為77×37×18,每個(gè)網(wǎng)格柵元的尺寸為2 m×2 m×2 m,涵蓋了待求解問題的全局范圍。圖2給出了在直接模擬情形下,計(jì)算時(shí)長分別為180 min、480 min和1 980 min時(shí)核電廠模型z=100 m平面Gamma通量分布的模擬結(jié)果。

從子圖2(a)可知,在直接模擬情形下,相當(dāng)一部分區(qū)域Gamma通量的計(jì)數(shù)結(jié)果為零,即出現(xiàn)了大量空網(wǎng)格的情形(圖中用白色表示),表明粒子很難甚至是無法輸運(yùn)到這部分深穿透區(qū)域。子圖2(b)和子圖2(c)是增加了模擬計(jì)算時(shí)長(或輸運(yùn)粒子歷史數(shù))后的結(jié)果,可以看出Gamma通量不為零的網(wǎng)格數(shù)有所增加,但計(jì)算時(shí)長達(dá)到1 980 min后,空網(wǎng)格仍然大量存在。因此,在直接蒙特卡羅模擬的情形下,若是簡單地通過增加模擬計(jì)算時(shí)長或輸運(yùn)粒子歷史數(shù)來改善空網(wǎng)格大量存在的不利局面,則其效果極其有限,從解決問題本身來說,產(chǎn)生的影響完全可以忽略不計(jì)。

圖2 Gamma通量分布的直接模擬結(jié)果Fig.2 The direct simulation results ofGamma flux distribution

圖3 基于蒙特卡羅正算迭代的全局減方差方法應(yīng)用流程Fig.3 The process of global variance reduction method based on Monte-Carlo forward transport iteration

為了改善甚至消除深穿透區(qū)域Gamma通量計(jì)數(shù)為零的現(xiàn)象(一般伴隨著相對誤差分布很不均勻的現(xiàn)象),進(jìn)而優(yōu)化蒙特卡羅模擬中計(jì)算資源在問題系統(tǒng)全局范圍內(nèi)的配置,本文利用基于蒙特卡羅正算輸運(yùn)迭代產(chǎn)生輸運(yùn)偏倚參數(shù)的方法,即公式(4),構(gòu)建了全局范圍內(nèi)的權(quán)重窗下限參數(shù)分布,進(jìn)而指導(dǎo)模擬粒子的輸運(yùn)過程,引導(dǎo)粒子盡可能地向問題系統(tǒng)的深穿透區(qū)域進(jìn)行輸運(yùn),通過增加深穿透區(qū)域的模擬粒子數(shù)和相應(yīng)地限制源區(qū)域及其附近的模擬粒子數(shù),來實(shí)現(xiàn)全局減方差和加速模擬結(jié)果收斂的目標(biāo),圖3給出了方法應(yīng)用的具體流程。圖4給出了基于網(wǎng)格通量的權(quán)重窗設(shè)置下,通過多次正算輸運(yùn)迭代得到的核電廠模型z=100 m平面Gamma通量分布的模擬結(jié)果和對應(yīng)的相對誤差(Re)分布。

圖4 Gamma通量迭代計(jì)算結(jié)果及相對誤差分布Fig.4 Gamma flux iteration calculation results and relative error distribution

由子圖4(b)可知,模型Gamma通量分布的模擬結(jié)果在經(jīng)過第一次迭代后就基本實(shí)現(xiàn)了消除深穿透區(qū)域零計(jì)數(shù)網(wǎng)格的目標(biāo)(即空網(wǎng)格率僅為0.04%),所耗費(fèi)的計(jì)算時(shí)長為480 min,與直接模擬時(shí)計(jì)算時(shí)長達(dá)到1 980 min仍存在大量零計(jì)數(shù)網(wǎng)格的情形相比,模擬質(zhì)量獲得了極大的提升,這一對比結(jié)果表明,利用網(wǎng)格通量構(gòu)建權(quán)重窗參數(shù)進(jìn)而指導(dǎo)下次迭代計(jì)算的粒子輸運(yùn)過程,有效增加了粒子輸運(yùn)到深穿透區(qū)域的機(jī)會,從計(jì)算資源在問題幾何各部分區(qū)域的分配方面來講,實(shí)現(xiàn)了計(jì)算資源的高效配置。在經(jīng)過第一次迭代后,深穿透區(qū)域網(wǎng)格計(jì)數(shù)結(jié)果的相對誤差與直接模擬相比,整體上得到了顯著的降低,且相對誤差為100%的網(wǎng)格數(shù)量也大量減少,表明蒙特卡羅模擬粒子在全局范圍分布不均勻的狀況相較直接模擬情形得到了顯著改善。由子圖4(d)可知,在經(jīng)過第三次迭代后,模型深穿透區(qū)域計(jì)數(shù)結(jié)果的相對誤差大幅度降低,即模擬結(jié)果的相對誤差分布在模型的全局范圍內(nèi)得到了展平,且絕大部分計(jì)數(shù)網(wǎng)格(約占總網(wǎng)格數(shù)的95%)的相對誤差都滿足了小于10%的統(tǒng)計(jì)要求,模擬質(zhì)量獲得了極大的提升,實(shí)現(xiàn)了全局減方差的目標(biāo)。經(jīng)過分析后發(fā)現(xiàn),相對誤差居高不下的網(wǎng)格(約占總網(wǎng)格數(shù)的5%)集中在密度較大的混凝土結(jié)構(gòu)中,這些結(jié)構(gòu)通常具備良好的屏蔽效果,通量水平較低,且引導(dǎo)的蒙特卡羅模擬粒子不易進(jìn)入其內(nèi)部,然而對于核電廠輻射環(huán)境分析而言,尤其針對人員可達(dá)性方面,需要的通量數(shù)據(jù)并不包括這些結(jié)構(gòu),因此產(chǎn)生的不利影響可以得到排除。

3 方法應(yīng)用分析

表1 總結(jié)了針對模型的蒙特卡羅模擬在迭代計(jì)算下獲得的減方差效果,表明相對于直接模擬,計(jì)算結(jié)果的FOM值均有所增加,最大增加到了直接情形的3倍。且其他參量如平均相對誤差的標(biāo)準(zhǔn)差σRe和空網(wǎng)格率也都朝著改善模擬質(zhì)量的方向發(fā)展。

表1 迭代計(jì)算的減方差效果數(shù)據(jù)指標(biāo)Table 1 The effect of global variance reduction technique based on the flux method

表2 中的數(shù)據(jù)表明,過低的權(quán)重窗上限參數(shù)的設(shè)置會限制一定時(shí)間內(nèi)源粒子抽樣數(shù)目,減少了系統(tǒng)內(nèi)實(shí)際模擬的粒子數(shù),從宏觀上來看,會增加網(wǎng)格計(jì)數(shù)結(jié)果的統(tǒng)計(jì)誤差,同時(shí),由于粒子發(fā)生過度分裂,導(dǎo)致低權(quán)重蒙特卡羅模擬粒子的輸運(yùn)占據(jù)了相當(dāng)一部分計(jì)算時(shí)間,影響了計(jì)算效率,使得權(quán)重窗技術(shù)的減方差效果沒有得到有效發(fā)揮,相關(guān)數(shù)據(jù)指標(biāo)如FOM等的計(jì)算值也反映了這一點(diǎn)。當(dāng)權(quán)重窗上限參數(shù)設(shè)置得過高時(shí)(如λ為1 000時(shí)),由于粒子很難發(fā)生分裂,導(dǎo)致系統(tǒng)內(nèi)高權(quán)重粒子仍然較多,使得網(wǎng)格計(jì)數(shù)結(jié)果的誤差并未得到有效降低。此外,當(dāng)λ為1 000時(shí),源粒子抽樣數(shù)目在相同的計(jì)算時(shí)間下和純模擬情形相同,說明權(quán)重窗技巧的應(yīng)用并未顯著增加系統(tǒng)內(nèi)蒙特卡羅模擬粒子的產(chǎn)生與輸運(yùn)機(jī)會,減方差的效果沒有得到有效實(shí)現(xiàn)。綜合表2中的數(shù)據(jù)指標(biāo)結(jié)果,對于本文針對的問題系統(tǒng),網(wǎng)格權(quán)重窗上限與下限參數(shù)之比設(shè)置為20～50時(shí)能夠有效地實(shí)現(xiàn)減方差技術(shù)的實(shí)際效果。

表2 不同上限參數(shù)設(shè)置下的減方差效果數(shù)據(jù)指標(biāo)(模擬時(shí)長為100 min)Table 2 The effect of global variance reduction technique under different upper weight bound(run time:100 min)

4 總結(jié)

本文基于核電廠嚴(yán)重事故下輻射環(huán)境的蒙特卡羅模擬需求,針對模擬計(jì)算中經(jīng)常存在的深穿透區(qū)域計(jì)數(shù)結(jié)果很難或者無法收斂的問題,闡述了全局減方差方法的思想及方法的應(yīng)用需要實(shí)現(xiàn)的目標(biāo),即展平蒙特卡羅模擬粒子在問題幾何的分布狀況和限制輸運(yùn)粒子的權(quán)重使其不出現(xiàn)過高和過低的狀況。在這兩個(gè)目標(biāo)的引導(dǎo)下,提出了構(gòu)建全局減方差偏倚參數(shù)的具體手段,即利用網(wǎng)格通量構(gòu)建問題幾何內(nèi)權(quán)重窗下限參數(shù)分布。在此理論基礎(chǔ)上,結(jié)合蒙特卡羅正算輸運(yùn)迭代的特點(diǎn)建立了迭代計(jì)算產(chǎn)生偏倚參數(shù)的具體方案,形成了一套完整的全局減方差方法的應(yīng)用流程。為了驗(yàn)證方法的應(yīng)用效果,本文針對核電廠模型開展了相關(guān)迭代計(jì)算,并與直接模擬的情形進(jìn)行了對比,對比結(jié)果表明,采用通量權(quán)窗指導(dǎo)粒子輸運(yùn),問題幾何的全局范圍內(nèi)都實(shí)現(xiàn)了減方差的效果,模擬質(zhì)量和效率(品質(zhì)因數(shù))都獲得了大幅度的提升,初步驗(yàn)證了方法在解決深穿透問題方面的有效性。最后,在對核電廠模型減方差的效果進(jìn)行具體量化的基礎(chǔ)上,進(jìn)一步分析了方法的應(yīng)用特點(diǎn)和實(shí)際效果,完善了對方法的實(shí)踐指導(dǎo)意義的研究。

致謝

本文的研究工作得到了國家科技重大專項(xiàng)

課題(2019ZX06004013,2017ZX06002003-001-002)的經(jīng)費(fèi)資助,特此感謝。