核電廠嚴重事故下輻射環境蒙特卡羅模擬及方法應用研究

王學棟,郭丁情,2,佟立麗,曹學武,*

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

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

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

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

1 相關概念和理論

1.1 蒙特卡羅方法的統計不確定性

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

式中:Sˉx——樣本均值ˉx的標準差;

S2——樣本方差;

N——模擬粒子歷史數。

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

品質因數的概念反映在計算效率上是其值與計算耗時之間的關系,可以作為一個表征參數,用來量化和比較蒙特卡羅模擬計算中減方差方法的實際應用效果。

1.2 全局減方差理論

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

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

式中:Wth,i——網格i內權重窗下限參數;

φi——網格i內的粒子通量計數結果;

λ——權重窗上限和下限參數之比,即Wsplit/Wth;

max(φ)——系統全局范圍內各個計數網格最大的粒子通量。

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

2 模擬計算與方法應用

2.1 計算模型

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

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

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

2.2 計算結果與討論

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

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

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

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

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

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

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

3 方法應用分析

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

表1 迭代計算的減方差效果數據指標Table 1 The effect of global variance reduction technique based on the flux method

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

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

4 總結

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

致謝

本文的研究工作得到了國家科技重大專項

課題(2019ZX06004013,2017ZX06002003-001-002)的經費資助,特此感謝。