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

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

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

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

1 方法描述

1.1 AIS方法

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

虛粒子位置為粒子沿輸運(yùn)方向與虛擬面的交點(diǎn)。權(quán)重w為碰撞粒子權(quán)重乘以碰撞粒子沿輸運(yùn)方向不經(jīng)歷碰撞直接抵達(dá)虛擬面的概率：

(1)

其中：w0為輸運(yùn)粒子權(quán)重；d為源或碰撞點(diǎn)到虛粒子位置的距離；Σt為總宏觀反應(yīng)截面。

圖1 AIS方法流程圖

1.2 MCShield程序

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

平行幾何在創(chuàng)建中可正常使用交、并、補(bǔ)及其他復(fù)雜幾何結(jié)構(gòu)，可與真實(shí)幾何重疊，因此理論上支持任意形狀的虛擬面，解決了原AIS方法虛擬面類(lèi)型受限的問(wèn)題，同時(shí)也增強(qiáng)了MCShield程序?qū)ζ帘螁?wèn)題的計(jì)算能力。

1.3 SDAIS方法

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

DXTRAN球上的虛粒子通過(guò)抽樣產(chǎn)生，并被賦予相應(yīng)的權(quán)重，以保證無(wú)偏性。圖3為DXTRAN球上的虛粒子抽樣方法。設(shè)粒子在P1發(fā)生碰撞，DXTRAN球心為P0，DXTRAN外球半徑為RO，內(nèi)球半徑為RI，P1P0長(zhǎng)度為L(zhǎng)。DXTRAN外球、內(nèi)球在P1P0方向上極角最大值分別為θO、θI，其余弦值分別為ηO、ηI。虛粒子位置PS在DXTRAN外球上，具體坐標(biāo)通過(guò)抽樣以P1P0方向?yàn)榛鶞?zhǔn)的極角和方位角確定。極角η余弦值按式(2)的概率密度函數(shù)p(η)抽樣。

圖2 SDAIS方法流程圖

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

(2)

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

虛粒子權(quán)重w為：

(3)

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

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

1) 對(duì)新產(chǎn)生的虛粒子，計(jì)算w/r2；

該算法不僅可應(yīng)用于SDAIS方法中，也可直接應(yīng)用于AIS方法，作為AIS方法統(tǒng)計(jì)小探測(cè)器的優(yōu)化選項(xiàng)。

2 算例驗(yàn)證

本文對(duì)美國(guó)NRC發(fā)布的NUREG/CR-6115 PWR注量率計(jì)算基準(zhǔn)題進(jìn)行了求解。NUREG/CR-6115 PWR基準(zhǔn)題由美國(guó)布魯克海文國(guó)家實(shí)驗(yàn)室(BNL)開(kāi)發(fā)。NUREG/CR-6115 PWR基準(zhǔn)題的反應(yīng)堆總熱功率為2 527.73 MW，堆芯含有204個(gè)燃耗深度不同的燃料組件。整個(gè)壓水堆結(jié)構(gòu)主要包括堆芯、吊籃、熱屏蔽層、壓力容器和混凝土生物屏蔽體等，其結(jié)構(gòu)如圖4所示?；炷辽锲帘误w范圍為R=335.915～549.275 cm，中子源強(qiáng)度為1.44×1020s-1 [18]。

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

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

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

采用FOM因子代表其計(jì)算效率，則：

(4)

其中：D為相對(duì)統(tǒng)計(jì)誤差；t為計(jì)算時(shí)間。FOM因子越大，代表計(jì)算效率越高。

3 結(jié)果分析

3.1 堆腔中子劑量?jī)x通量(E>0.1 MeV)及能譜

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

統(tǒng)計(jì)能量0.1 MeV以上的中子通量，結(jié)果列于表1。由表1可見(jiàn)，以DORT程序的結(jié)果作為參考基準(zhǔn)，幾種方法結(jié)果均符合得較好，與DORT程序的相對(duì)偏差均在5%以?xún)?nèi)。AIS方法的計(jì)算效率是IMP-MC的8倍，而本文提出的SDAIS方法的計(jì)算效率約是IMP-MC方法的56倍，相比AIS方法提高了7倍。

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

表1 堆腔中子劑量?jī)x通量

圖5 堆腔劑量?jī)x中子能譜

圖6 堆腔劑量?jī)x中子能譜FOM因子

3.2 混凝土生物屏蔽體內(nèi)球形探測(cè)器通量

在12個(gè)探測(cè)器、5種計(jì)算方法的通量結(jié)果中，除直徑1 cm探測(cè)器的AIS方法和AIS*方法的統(tǒng)計(jì)誤差較大外，其他結(jié)果的統(tǒng)計(jì)誤差均在10%以下。直徑5 cm的6個(gè)探測(cè)器的歸一化通量結(jié)果如圖7所示，各方法的計(jì)算結(jié)果符合得較好，各方法間的平均偏差在10%以?xún)?nèi)。

圖7 生物屏蔽體內(nèi)5 cm球形探測(cè)器中子通量

圖8 生物屏蔽體內(nèi)球形探測(cè)器FOM因子比較

對(duì)不同位置探測(cè)器的FOM因子進(jìn)行平均，結(jié)果如圖8所示。對(duì)于探測(cè)器直徑5 cm的情形，AIS*方法效率最高，SDAIS方法次之；對(duì)于探測(cè)器直徑1 cm的情形，SDAIS方法效率最高，傳統(tǒng)AIS方法最差。這是因?yàn)樘綔y(cè)器的直徑與直接進(jìn)入探測(cè)器區(qū)域的粒子數(shù)呈正相關(guān)，而后者直接決定了統(tǒng)計(jì)誤差。當(dāng)探測(cè)器直徑較小時(shí)，直接進(jìn)入探測(cè)器區(qū)域的粒子數(shù)很少，此時(shí)需使用指向概率的方法進(jìn)行統(tǒng)計(jì)。此外可看到，使用基于探測(cè)器位置的虛粒子賭分裂算法的計(jì)算效率比傳統(tǒng)方法高約1個(gè)量級(jí)。另外，SDAIS方法和F5AIS*方法均使用了指向概率方法，但SDAIS方法比F5AIS*方法的計(jì)算效率高些，這是因?yàn)镾DAIS方法僅在最后1個(gè)虛擬面外才進(jìn)行指向概率的計(jì)算，計(jì)算量較小。

4 結(jié)論

本文在AIS方法的基礎(chǔ)上，針對(duì)小探測(cè)器問(wèn)題，優(yōu)化了虛粒子賭分裂算法，將AIS方法和指向概率方法相結(jié)合，提出了SDAIS方法，在自主開(kāi)發(fā)的蒙特卡羅程序MCShield中進(jìn)行了實(shí)現(xiàn)。使用NUREG/CR-6115 PWR基準(zhǔn)題對(duì)該方法的正確性及計(jì)算效率進(jìn)行了驗(yàn)證。計(jì)算結(jié)果顯示，當(dāng)探測(cè)器較小時(shí)，SDAIS方法的計(jì)算效率比傳統(tǒng)AIS方法高約1個(gè)量級(jí)，比常規(guī)蒙特卡羅方法高約2個(gè)量級(jí)?；谔綔y(cè)器位置的虛粒子賭分裂算法可較為有效地提升傳統(tǒng)AIS方法在小探測(cè)器問(wèn)題中的計(jì)算效率。