基于自動重要抽樣方法的減方差技巧體系構建與驗證

武 禎,郝以昇,浦彥恒,周 揚,杲申申,邱 睿,*,馬銳垚,張 輝,李君利

(1.清華大學 工程物理系,北京 100084;2.粒子技術與輻射成像教育部重點實驗室,北京 100084;3.同方威視技術股份有限公司,北京 100084)

輻射屏蔽系統的設計與分析對核設施的安全具有重要意義,隨著核電技術的發展,對精細化屏蔽計算的要求越來越高。使用蒙特卡羅方法進行屏蔽計算可準確地模擬粒子的行為,與其他方法相比,其在計算精度、物理模型的完備性、幾何描述能力等方面具有無可比擬的優勢[1]。但蒙特卡羅方法也存在收斂速度慢、計算較費時的缺點,這成為其在實際工程中應用的主要障礙。Wagner等[1]按照計算目標將粒子輸運問題分為源-探測器問題、區域問題和全局問題,為解決這3類問題中存在的深穿透問題,研究人員提出了一系列減方差技巧。按照實現原理,目前已有的減方差技巧可分為偏倚、半解析和耦合減方差方法[2]。幾何重要性方法和權窗方法是使用廣泛的偏倚減方差方法,它們通過人為修改粒子權重和數量減小方差。半解析減方差方法包括隱式俘獲法、統計估計法、點通量法和DXTRAN球法等,它們主要通過在粒子輸運中運用解析方法降低方差,然而其并沒有改變粒子的隨機游動過程,因此對于深穿透問題效率提升不大。耦合減方差方法的特點是進行多次計算,其中一類是在空間上進行耦合,例如確定論和蒙特卡羅的空間耦合[3-4],其利用確定論方法計算大屏蔽、幾何簡單部分,而利用蒙特卡羅方法計算源項、幾何復雜的部分;又例如蒙特卡羅空間分層面源續算方法[5-8]也屬于空間耦合方法。另一類是利用一次計算給出合理的幾何重要性或權窗參數,再利用幾何重要性或權窗方法進行二次計算。根據一次計算的方法可分為蒙特卡羅正算、蒙特卡羅伴隨、確定論正算和確定論伴隨。一次計算采用正算的耦合減方差方法適用于全局問題和部分區域問題,利用正算通量與權重的正比關系生成權窗參數。文獻[8]用確定論程序進行一次計算正算生成空間上的權窗參數,文獻[9]將其擴展為空間-能量網格。采用蒙特卡羅程序進行一次計算正算的方法有MAGIC方法[10]和基于蒙特卡羅正算的全局減方差方法[11]等。一次計算采用確定論伴隨計算的一致性共軛驅動重要抽樣(CADIS)[12]方法最初適用于源-探測器問題,改進后的前向加權一致性伴隨驅動重要抽樣(FW-CADIS)[1]方法適用于區域問題和全局問題。一次計算采用蒙特卡羅伴隨計算的耦合減方差方法方面,邱有恒等[13]曾仿照CADIS方法在MCNP上使用多群伴隨蒙特卡羅計算生成源偏倚和權窗參數。

為了解決深穿透問題[14-15],清華大學輻射防護與環境保護實驗室提出了自動重要抽樣(AIS)方法的減方差技巧[16-19],且在為輻射屏蔽計算設計開發的中子/光子/電子耦合輸運蒙特卡羅程序MCShield[20]中進行了實現。AIS方法引入基于重要性采樣和統計估計方法的虛擬面,將空間劃分為多個子空間層,虛粒子在虛擬面上產生并傳送到下一個子空間,隨后,進行自動粒子權重調整和數量控制,使虛擬面上的虛粒子數與源粒子數相等。AIS方法需要在源項與統計區域之間定義一系列虛擬面。在MCShield程序中,用戶可根據實際需求,沿特定方向設置平面、圓柱面和球面3種類型的虛擬面,同時需要保證虛擬面在計算區域內不相交。雖然AIS方法在解決區域問題時取得了良好的效果[2],但在源-探測器問題和全局問題應用時效率還不是很高,因此MCShield研究團隊后續又基于AIS方法提出了一系列減方差技巧用于解決這兩類問題。本文以此為基礎構建基于AIS方法的減方差技巧體系,并對其中的部分減方差技巧進行驗證。

1 基于AIS方法的減方差技巧體系

屏蔽計算問題根據計算目標不同可分為源-探測器問題、區域問題和全局問題。針對不同的目標,適用的減方差方法不同,MCShield研究團隊針對這3類問題提出了基于AIS方法的系列減方差技巧。源-探測器問題通常分為小體積探測器問題和迷宮孔道兩類問題,針對小體積探測器問題,提出了小探測器自動重要抽樣(SDAIS)方法和AIS伴隨蒙特卡羅的耦合減方差(AIS-CADIS)方法。SDAIS方法基于探測器位置進行虛粒子數量控制并將AIS方法與指向概率結合,AIS-CADIS方法將AIS方法引入伴隨蒙特卡羅輸運計算中,利用AIS伴隨生成源偏倚和權窗參數減小探測器的統計誤差[20-21]。針對迷宮孔道問題,提出了FPAIS方法[2,22],采用方向概率引導粒子輸運,以實現減方差的目的。

在區域問題中,往往需要處理大體積和多個探測器。針對區域問題,提出了CNP-AIS方法[2],采用統計估計法和分層輸運法實現中子-光子耦合輸運。針對全局問題,涉及較大的空間尺度,提出了網格化-AIS方法[20],使整個空間的粒子密度盡可能均勻,以獲得統一的全局統計誤差。以上述系列減方差技巧為基礎,本文構建了基于AIS方法的減方差技巧體系,如圖1所示。

2 源-探測器問題減方差技巧驗證

2.1 小探測器自動重要抽樣方法驗證

MCShield研究團隊采用NUREG/CR-6115 PWR壓力容器計算基準題[23]對SDAIS方法進行驗證,選取該基準題的標準核燃料組件布局方案,計算反應堆堆腔內中子劑量儀的中子通量。反應堆中子劑量計位于反應堆腔體內徑向距離r=319.56～320.56 cm和軸向高度z=176.27～178.27 cm處,如圖2所示。

圖2 NUREG/CR-6115 PWR基準題幾何示意圖Fig.2 NUREG/CR-6115 PWR benchmark geometry diagram

分別采用MCShield程序的幾何重要性法(IMP-MC)、AIS方法和SDAIS方法進行計算,比較3種方法的正確性和效率。其中,AIS方法使用4個圓柱形中子虛擬面,半徑分別為188、215、230和300 cm。SDAIS方法與AIS方法相比,除了使用半徑分別為188、215、230 cm的3個圓柱中子虛擬面外,還進行兩點改進:1) 考慮虛粒子位置對探測器貢獻的影響,添加修正因子1/d2(d為虛粒子距離探測器的距離)對虛粒子的數量進行控制,達到將更多的粒子向探測器區域輸運的目的;2) 在探測器所在位置處設置一個DXTRAN球,中心坐標為(317.04 cm,43.85 cm,177.27 cm),外球半徑為7 cm,內球半徑為6.5 cm,增加粒子到達探測器區域附近并在其中發生碰撞的概率。

本算例中,選取基準題中確定論程序DORT作為基準參考結果,使用品質因子(FOM)來評估計算效率,FOM定義為:

FOM=1/(R2T)

(1)

其中:R為統計誤差;T為計算時間。

堆腔中子劑量儀通量(E>0.1 MeV)列于表1,其中,歸一化通量表示由每個粒子引起的通量。可看出,中子通量在中子劑量儀中的衰減約為10-10cm-2。上述方法的計算結果與DORT程序結果基本一致,SDAIS方法的計算效率約為IMP-MC的56倍、AIS方法的7倍。事實證明,SDAIS方法對于源-探測器的問題具有更高的計算效率。除對反應堆堆腔的中子劑量儀的中子通量進行驗證外,還驗證了混凝土屏蔽層內的中子通量[22]。

表1 堆腔中子劑量儀通量Table 1 Reactor cavity neutron dosimeter flux

2.2 AIS伴隨蒙特卡羅方法驗證

基于伴隨計算的耦合減方差方法在源-探測器問題中具有良好的效果。目前基于伴隨計算的耦合減方差方法中,主流的方法是在一次計算中使用確定性伴隨計算的CADIS方法。基于蒙特卡羅方法雖也可進行伴隨計算,但在深穿透問題中基于蒙特卡羅方法正向計算的伴隨問題仍是深穿透問題。為解決這一問題,本文提出了一種基于AIS伴隨蒙特卡羅方法的耦合減方差方法,將AIS方法引入伴隨蒙特卡羅方法中,以解決深穿透問題中收斂速度慢的問題。

本文將基于AIS伴隨蒙特卡羅的耦合減方差方法應用到某商用反應堆屏蔽算例的計算中。某商用反應堆屏蔽算例的二維模型如圖3所示,二維切面圖的切面均經過測點。探測器位于主管道的保溫層內部,共有4個測點,編號為P1、P2、P3、P4,統計各測點快中子(E>1 MeV)和熱中子(E<0.6 eV)通量。使用MCNP原始伴隨程序和AIS伴隨程序分別生成源偏倚和權窗參數,CADIS方法中確定論方法計算的源偏倚和權窗參數作為參考。在AIS伴隨方法中采取一系列球面作為虛擬面,如以P1點作為源項的伴隨計算時,AIS虛擬面設置為以坐標(0 cm,-40 cm,30 cm)為球心的一系列球面,半徑依次為380、320、260、200、140、100、40 cm。

圖3 某商用反應堆屏蔽算例的二維模型Fig.3 Two-dimensional model of commercial reactor shielding example

以測量值為基準,統計了3種方法的計算結果與測量值的相對偏差,如圖4所示。其中,AIS-CADIS代表基于AIS伴隨蒙特卡羅計算的正算結果,MC-CADIS代表基于原始伴隨蒙特卡羅計算的正算結果,CADIS代表使用CADIS方法的JMCT程序[24]計算結果。可看出,CADIS和AIS-CADIS在快中子通量結果上較為接近,最大相對偏差小于5%,與測量值的相對偏差也均在20%以內;在熱中子通量結果上,CADIS和AIS-CADIS最大相對偏差約為20%;MC-CADIS結果較CADIS和AIS-CADIS的結果偏小超過20%,遠大于相對統計誤差,原因是伴隨通量誤差較大,導致計算偏倚源項時出現系統誤差。這說明基于原始伴隨蒙特卡羅的減方差方法在深穿透問題中有一定的局限性。

圖4 快中子(a)和熱中子(b)通量與測量值相對偏差Fig.4 Relative deviation of fast neutron (a) and thermal neutron (b) fluxes from measured value

圖5為快中子和熱中子通量FOM的比較。可看出,對于快中子通量結果,CADIS的計算效率約為AIS-CADIS的1～3倍。對于熱中子通量結果,AIS-CADIS的計算效率約為CADIS的5～32倍。總體上說,AIS-CADIS的計算效率略優于CADIS。

圖5 快中子(a)和熱中子(b)通量FOM比較Fig.5 Comparison of fast neutron (a) and thermal neutron (b) flux FOM

3 全局問題減方差技巧驗證

本文通過簡化的反應堆屏蔽計算算例驗證了網格化-AIS方法,圖6為自設反應堆屏蔽算例幾何示意圖。

圖6 簡化反應堆屏蔽算例幾何示意圖Fig.6 Geometry diagram of simplified reactor shielding calculation example

反應堆堆芯是一半徑為0.5 m、高為1 m的圓柱體,其下表面與水層底部共面。防護罩外側為半徑為1.5 m、高為4 m的水層,外側為真空邊界的1 m混凝土層。兩根0.1 m半徑的空氣管分別在徑向和軸向穿過水層和混凝土層。中子源各向同性,能量分布符合瓦特裂變譜。本算例測試了直接蒙特卡羅方法、AIS方法和網格化-AIS方法在各向異性較強幾何下的計算表現,3種方法計算得到的中子通量及相對統計誤差三維分布沿空氣管道中心處的切面如圖7～9所示。表2列出了簡化反應堆屏蔽算例計算結果。網格化-AIS方法的平均相對統計誤差為3.3%,在實際屏蔽計算中,通常要求相對統計誤差需要小于5%。另外,網格化-AIS方法的FOM是AIS方法的12倍左右,是直接蒙特卡羅方法的290倍左右。

表2 簡化反應堆屏蔽算例計算結果Table 2 Calculation result of simplified reactor shielding example

圖7 直接蒙特卡羅方法的通量(a)及相對統計誤差(b)Fig.7 Analog Monte Carlo method flux (a) and relative statistical error (b)

圖8 AIS方法的通量(a)及相對統計誤差(b)Fig.8 AIS method flux (a) and relative statistical error (b)

圖9 網格化-AIS方法的通量(a)及相對統計誤差(b)Fig.9 Grid-AIS method flux (a) and relative statistical error (b)

4 區域問題減方差技巧改進

對于區域問題使用的減方差方法CNP-AIS,需要在源項與統計區域之間定義一系列的平面、圓柱面或球面虛擬面,將整個計算區域劃分為多個子空間。在實際應用中,虛擬面需要用戶憑借經驗進行設置,計算結果的準確性和效率取決于虛擬面設置的好壞。在不同的幾何條件下,虛擬面的設置方法及最優位置的選擇是不同的,這也就要求用戶具備較為豐富的屏蔽計算經驗,不利于程序的通用性。同時,目前AIS方法僅支持平面、圓柱面或球面3種虛擬面,不足以適應各種復雜幾何的情況。

為解決以上問題,本文結合CADIS方法,提出了適用于復雜幾何的虛擬面自動生成與調整算法,其基本過程為:1) 讀取各CADIS方法計算的相空間網格參數,得到計算空間中各網格的響應貢獻分布情況;2) 基于網格的響應貢獻分布生成源偏倚參數,同時繪制響應貢獻分布等值面,基于等值面自動生成多個網格虛擬面,用于后續的AIS方法虛粒子生成;3) 將虛擬面所在網格中的能量權重結果進行匯總,生成該虛擬面的能量權重參數,用于對虛擬面上的虛粒子進行能量權重修正;4) 在對到達虛擬面的虛粒子進行再抽樣時,需要將該虛粒子自身的權重乘以其對應能量權重參數,再對其進行輪盤賭或分裂生成新的虛粒子,然后作為下一個子空間的源粒子進行繼續輸運。

5 結論

本文對MCShield研究團隊針對屏蔽計算深穿透問題研發的基于AIS方法的系列減方差技巧進行了總結,構建了基于AIS方法的減方差技巧體系,用于解決區域問題、源-探測器問題和全局問題,并對其中的一些方法進行了介紹及驗證。采用NUREG/CR-6115 PWR壓力容器計算基準題驗證了SDAIS方法,結果表明SDAIS方法對于源-探測器問題的計算效率高于AIS方法。提出了AIS-CADIS方法,將AIS方法引入到蒙特卡羅伴隨計算中,對某商用反應堆主管道內的測點通量進行了計算驗證,結果表明,該方法與CADIS方法的效率相當,且僅使用一套程序,避免了確定論與蒙特卡羅耦合計算帶來的問題。針對全局問題,使用簡化的反應堆屏蔽計算算例驗證了網格化-AIS方法,結果表明,與AIS方法和直接蒙特卡羅方法相比,網格化-AIS方法可有效提升計算效率。