用于厚屏蔽小探測器的蒙特卡羅模擬減方差方法研究

杲申申，李君利，武 禎，馬銳垚，王 鑫，邱 睿，李春艷，張 輝

(1.清華大學 工程物理系，北京 100084；2.粒子技術與輻射成像教育部重點實驗室，北京 100084；3.同方威視技術股份有限公司，北京 100084；4.北京應用物理與計算數學研究所，北京 100094)

蒙特卡羅方法是屏蔽計算的重要工具，它可精確描述三維復雜幾何與物理過程，較為容易地給出計算誤差。但蒙特卡羅方法收斂速度較慢，在厚屏蔽問題及小探測器問題中統計誤差較大。厚屏蔽問題的特點是源和統計區域間的屏蔽體足夠厚，以至于粒子穿透屏蔽體到達統計區域的概率極小，使用常規蒙特卡羅輸運統計誤差較大。在小探測器問題中，探測器體積小導致粒子輸運到探測器中的概率小，同樣會增大蒙特卡羅統計誤差。在反應堆屏蔽計算中，厚屏蔽和小探測器問題經常同時存在，常規蒙特卡羅方法難以有效解決。近年來，隨著計算機性能的提升和高性能計算科學的發展，大規模并行計算可模擬海量的粒子數，在一定程度上緩解了深穿透問題帶來的計算壓力。同時，國內外學者在減方差技巧方面也做了大量的研究工作。在實際應用中，常用的減方差技巧有幾何重要性、網格權窗、指數變換及針對小探測器的點探測器、DXTRAN方法等[1-2]，其中，確定論伴隨-蒙特卡羅耦合計算方式是目前屏蔽計算減方差的主要方法之一。該方法的流程是先由確定論方法求解原問題的伴隨方程，生成蒙特卡羅計算的網格權窗參數，再用蒙特卡羅程序進行計算。但該方法需額外重建幾何并進行兩次輸運計算，且確定論程序對復雜幾何的描述也不夠精確[3-6]。

2004年，文獻[7]針對厚屏蔽、深穿透問題提出了自動重要抽樣(AIS)方法；文獻[8-9]運用該方法在內照射劑量計算、探測器校正因子計算等方面獲得了較好的效果；文獻[10-14]從效率和實用性等方面對該方法進行了改進。AIS方法基于重要抽樣、統計估計法和分層輸運的思想，無需設置重要性或權窗參數，正確性和計算效率在一系列基準題上得到了驗證，效率相比直接模擬可提高1～2個數量級。對厚屏蔽體外小體積探測器問題，AIS方法雖可生成虛粒子，使較多的粒子到達屏蔽體外，但因為探測器體積較小，實際抵達探測器的粒子仍然不足。而基于指向概率的DXTRAN方法[15]對小探測器的估計更為準確。因此本文基于AIS方法和指向概率方法，提出并實現了小探測器自動重要抽樣(SDAIS)方法，同時針對小探測器問題，對AIS方法原有的虛粒子賭分裂算法進行改進，提出基于探測器位置的虛粒子賭分裂算法，使用NUREG/CR-6115 PWR基準題算例，對該方法應用于厚屏蔽小探測器問題的能力進行驗證。

1 方法描述

1.1 AIS方法

AIS方法自2004年提出以來，經多次改進，在反應堆屏蔽、人體體素模型計算等問題的應用中取得了較好的效果。該方法基于重要性抽樣和統計估計法，引入虛擬面將空間劃分為多層子空間，在虛擬面上產生虛粒子向下一層子空間輸運，并進行自動的粒子權重調整和數量控制。圖1為AIS方法流程圖，其流程為：1) 引入K個虛擬面，輸運空間劃分為K+1個子空間；2) 當源項抽樣或粒子發生碰撞時，在虛擬面上生成虛粒子，并維持虛擬面上的虛粒子數等于源粒子數；3) 粒子若在輸運中到達當前虛擬面，則將其殺死；4) 當前虛擬面上的虛粒子作為下一子空間的源粒子繼續輸運。

虛粒子位置為粒子沿輸運方向與虛擬面的交點。權重w為碰撞粒子權重乘以碰撞粒子沿輸運方向不經歷碰撞直接抵達虛擬面的概率：

(1)

其中：w0為輸運粒子權重；d為源或碰撞點到虛粒子位置的距離；Σt為總宏觀反應截面。

圖1 AIS方法流程圖

1.2 MCShield程序

AIS方法由范佳錦[7]在MCNP上進行了實現，虛擬面的設置使用RDUM卡，與普通幾何設置獨立，需額外編寫代碼判斷粒子是否會穿過虛擬面，因此虛擬面的類型受限，僅支持平面、球面、圓柱面等簡單幾何。MCShield程序[13,16]是清華大學工程物理系輻射防護實驗室自主開發的中子、光子、電子耦合輸運輻射屏蔽蒙特卡羅計算程序，采用CSG幾何架構，支持CAD模型幾何導入和計算結果三維可視化，具有強大的前、后處理能力。本文將AIS方法在MCShield程序中進行實現和改進，虛擬面的設置利用Geant4[17]的平行幾何思想進行實現。平行幾何是獨立于真實幾何的一套幾何系統，只有幾何形狀信息，沒有材料等信息。AIS方法需在粒子每次碰撞時判斷與虛擬面的距離，且虛粒子需在虛擬面駐留。在平行幾何中，粒子在真實幾何中進行幾何邊界檢查、駐留及碰撞，在平行幾何中僅做幾何邊界檢查和駐留。因此，平行幾何能滿足算法實現的需求，又可盡量減少對常規輸運的效率影響。

平行幾何在創建中可正常使用交、并、補及其他復雜幾何結構，可與真實幾何重疊，因此理論上支持任意形狀的虛擬面，解決了原AIS方法虛擬面類型受限的問題，同時也增強了MCShield程序對屏蔽問題的計算能力。

1.3 SDAIS方法

對于厚屏蔽體外小探測器問題，AIS方法未利用探測器信息，產生的虛粒子并不會向探測器聚集，而是大致在虛擬面上均勻分布。當探測器體積較小時，實際抵達探測器的粒子仍然不足。在這種情況下，DXTRAN方法對小探測器的估計更為準確。因此本文將AIS方法和DXTRAN方法結合，提出SDAIS方法。該方法在源與探測器之間建立K個AIS虛擬面，將空間劃分為K+1個子空間，探測器位于第K+1個，即最后1個子空間中。在最后1個子空間中，設置DXTRAN球覆蓋探測器。粒子在前K層子空間輸運時遵循AIS方法，在最后1個子空間中輸運時，每當粒子發生碰撞，就會在DXTRAN球上產生虛粒子，使更多的粒子輸運到探測器周圍。圖2為SDAIS方法的流程圖，其具體流程為：1) 當最后1層AIS虛擬面外的粒子發生碰撞時，在DXTRAN球上產生虛粒子；2) 若普通粒子進入DXTRAN球則將其殺死；3) 當普通粒子輸運完畢后，再輸運DXTRAN球上的虛粒子。若粒子穿出DXTRAN球則進行輪盤賭，存活下來的粒子轉化為普通粒子繼續輸運。

DXTRAN球上的虛粒子通過抽樣產生，并被賦予相應的權重，以保證無偏性。圖3為DXTRAN球上的虛粒子抽樣方法。設粒子在P1發生碰撞，DXTRAN球心為P0，DXTRAN外球半徑為RO，內球半徑為RI，P1P0長度為L。DXTRAN外球、內球在P1P0方向上極角最大值分別為θO、θI，其余弦值分別為ηO、ηI。虛粒子位置PS在DXTRAN外球上，具體坐標通過抽樣以P1P0方向為基準的極角和方位角確定。極角η余弦值按式(2)的概率密度函數p(η)抽樣。

圖2 SDAIS方法流程圖

圖3 DXTRAN球上虛粒子的抽樣方法

(2)

其中：Q為調整系數，在該算法中取5；U=Q(1-ηI)+ηI-ηO。方位角按各向同性抽樣。

虛粒子權重w為：

(3)

其中，p(μ)為粒子碰撞散射角余弦為μ的概率密度。

1.4 基于探測器位置的虛粒子賭分裂算法

1) 對新產生的虛粒子，計算w/r2；

該算法不僅可應用于SDAIS方法中，也可直接應用于AIS方法，作為AIS方法統計小探測器的優化選項。

2 算例驗證

本文對美國NRC發布的NUREG/CR-6115 PWR注量率計算基準題進行了求解。NUREG/CR-6115 PWR基準題由美國布魯克海文國家實驗室(BNL)開發。NUREG/CR-6115 PWR基準題的反應堆總熱功率為2 527.73 MW，堆芯含有204個燃耗深度不同的燃料組件。整個壓水堆結構主要包括堆芯、吊籃、熱屏蔽層、壓力容器和混凝土生物屏蔽體等，其結構如圖4所示。混凝土生物屏蔽體范圍為R=335.915～549.275 cm，中子源強度為1.44×1020s-1 [18]。

a——全堆模型俯視圖；b——θ=0°處的切面視圖

本文選取基準題中的標準核燃料組件布置方案作為計算對象，計算了參考文獻[18]中的堆腔中子劑量儀能譜，堆腔中子劑量儀位于堆腔內部R=319.56～320.56 cm、Z=176.27～178.27 cm處(圖4)。為了比較計算的正確性與計算效率，分別應用MCShield幾何分裂和輪盤賭(IMP-MC)方法、AIS方法及SDAIS方法進行了計算，確定論SN程序DORT的計算結果作為參考解。

同時為進一步驗證SDAIS方法在厚屏蔽情況下的計算效率，本文在距堆芯更遠的混凝土內部R=370、400、430、460、490、520 cm，Z=191.06 cm處各設置了直徑為1 cm和5 cm的球形探測器，共12個。使用AIS方法、SDAIS方法及點探測器統計的AIS方法(F5AIS)進行了計算。同時為驗證基于探測器位置的虛粒子賭分裂算法的通用性，也使用基于該算法的AIS方法和點探測器統計的AIS方法(分別記為AIS*和F5AIS*)進行了計算。

采用FOM因子代表其計算效率，則：

(4)

其中：D為相對統計誤差；t為計算時間。FOM因子越大，代表計算效率越高。

3 結果分析

3.1 堆腔中子劑量儀通量(E>0.1 MeV)及能譜

IMP-MC方法的樣本數為4×108，計算時間為13 736 min。AIS方法引入了4個圓柱形中子虛擬面，半徑分別為188、215、230、300 cm，樣本數為5×106，計算時間為206 min。SDAIS方法引入了3個半徑分別為188、215、230 cm的圓柱形中子虛擬面和1個DXTRAN球。DXTRAN球的中心為劑量儀中心(317.04 cm，43.85 cm，177.27 cm)，外球半徑7 cm，內球半徑6.5 cm，使用了基于探測器位置的虛粒子賭分裂算法，樣本數為5×106，計算時間為251 min。

統計能量0.1 MeV以上的中子通量，結果列于表1。由表1可見，以DORT程序的結果作為參考基準，幾種方法結果均符合得較好，與DORT程序的相對偏差均在5%以內。AIS方法的計算效率是IMP-MC的8倍，而本文提出的SDAIS方法的計算效率約是IMP-MC方法的56倍，相比AIS方法提高了7倍。

采用與基準題中DORT方法一致的47個能量間隔，統計堆腔中子劑量儀的中子能譜。IMP-MC方法、AIS方法和SDAIS方法的中子能譜計算結果及FOM因子如圖5、6所示。IMP-MC方法的中子能譜中有統計值的平均蒙特卡羅相對統計誤差為26.0%，平均FOM因子為0.001 08 min-1；AIS方法的中子能譜中有統計值的平均蒙特卡羅相對統計誤差是32.9%，平均FOM因子為0.044 8 min-1；SDAIS方法的中子能譜中有統計值的平均蒙特卡羅相對統計誤差為7.66%，平均FOM因子為0.679 min-1。除AIS方法因統計誤差較大導致與其他計算結果偏離較大，其他計算結果符合得較好。SDAIS方法的平均計算效率比AIS方法高約1個量級，比IMP-MC方法高約2個兩級。

表1 堆腔中子劑量儀通量

圖5 堆腔劑量儀中子能譜

圖6 堆腔劑量儀中子能譜FOM因子

3.2 混凝土生物屏蔽體內球形探測器通量

在12個探測器、5種計算方法的通量結果中，除直徑1 cm探測器的AIS方法和AIS*方法的統計誤差較大外，其他結果的統計誤差均在10%以下。直徑5 cm的6個探測器的歸一化通量結果如圖7所示，各方法的計算結果符合得較好，各方法間的平均偏差在10%以內。

圖7 生物屏蔽體內5 cm球形探測器中子通量

圖8 生物屏蔽體內球形探測器FOM因子比較

對不同位置探測器的FOM因子進行平均，結果如圖8所示。對于探測器直徑5 cm的情形，AIS*方法效率最高，SDAIS方法次之；對于探測器直徑1 cm的情形，SDAIS方法效率最高，傳統AIS方法最差。這是因為探測器的直徑與直接進入探測器區域的粒子數呈正相關，而后者直接決定了統計誤差。當探測器直徑較小時，直接進入探測器區域的粒子數很少，此時需使用指向概率的方法進行統計。此外可看到，使用基于探測器位置的虛粒子賭分裂算法的計算效率比傳統方法高約1個量級。另外，SDAIS方法和F5AIS*方法均使用了指向概率方法，但SDAIS方法比F5AIS*方法的計算效率高些，這是因為SDAIS方法僅在最后1個虛擬面外才進行指向概率的計算，計算量較小。

4 結論

本文在AIS方法的基礎上，針對小探測器問題，優化了虛粒子賭分裂算法，將AIS方法和指向概率方法相結合，提出了SDAIS方法，在自主開發的蒙特卡羅程序MCShield中進行了實現。使用NUREG/CR-6115 PWR基準題對該方法的正確性及計算效率進行了驗證。計算結果顯示，當探測器較小時，SDAIS方法的計算效率比傳統AIS方法高約1個量級，比常規蒙特卡羅方法高約2個量級。基于探測器位置的虛粒子賭分裂算法可較為有效地提升傳統AIS方法在小探測器問題中的計算效率。