300＃研究堆安全棒中子注量率計算中的減方差方法對比及應用

唐鳳平，劉耀光，楊萬奎，楊 鑫，2

（1.中國工程物理研究院 核物理與化學研究所，四川 綿陽 621900；2.清華大學 工程物理系，北京 100084）

300＃研究堆的安全棒為碳化硼吸收體，為研究碳化硼在整個反應堆壽期內的輻照性能，需要獲取整根碳化硼吸收體沿高度方向的中子注量率分布情況。但安全棒頂端離堆芯較遠，且有效吸收體體積較小，此類問題即屬于遠距離小體積計數，經計算，直接使用徑跡長度或面通量計數根本得不到安全棒頂端的計算結果，因此需要考慮減方差方法［1］。類似的小體積遠距離輸運問題是蒙特卡羅方法應用中的一個經典問題，這是由粒子到達計數區域附近的概率較小所致。清華大學已開展內照射小器官劑量計算中減方差技巧的比較和應用［2］，通過此方法獲得了較好的結果。本文針對該安全棒頂端的計數問題開展若干減方差方法的對比研究。

1 堆芯布置及安全棒位置

300＃研究堆額定功率3 MW，是一座游泳池式研究堆，堆芯包括活性區（燃料組件）、反射層（鈹組件和石墨組件）以及垂直孔道。堆芯布置如圖1所示［3］。堆芯中心為燃料組件，鈹組件和石墨組件作為反射層。其中燃料組件22個（16根燃料元件的燃料組件12個，15根燃料元件的燃料組件10個），鈹組件30個，石墨組件33個，控制棒12根（位于石墨盒的控制棒導管內和15根燃料元件的元件盒內），垂直空腔孔道1個，垂直水腔孔道4個。

圖1 反應堆堆芯布置及安全棒位置示意圖Fig.1 Schematic diagram of reactor core loading and safety rod's location

安全棒為碳化硼吸收體，有效吸收尺寸為φ17.5mm×600 mm。1AB 代表1號安全棒，2AB代表2號安全棒，1AB位于中央孔道的左側，2AB位于中央孔道的右側，其高度位置示意圖如圖2所示，即反應堆啟動時將兩根安全棒提升至反應堆頂部，之后再提升其余棒使反應堆達臨界。

圖2 安全棒高度位置示意圖Fig.2 Schematic diagram of safety rod's location in height

2 減方差原理

2.1 一般原理

按減方差原理，可將若干減方差方法分為4類［1］：

1）截斷方法

該類方法是最簡單直接的方法，主要包括幾何截斷、能量截斷和時間截斷。

2）數量控制方法

該類方法通過粒子分裂和俄羅斯輪盤賭來控制抽樣數量。在重要的區域增加抽樣，但每個粒子的權重降低；在不重要的區域減少抽樣，但每個粒子的權重增大。該方法通過調整權重以保證結果無偏。主要包括幾何、能量、時間相關的賭分裂、權重截斷以及權窗。

3）修正抽樣方法

該類方法通過改變統計抽樣以增加每個粒子的計數率。對任意隨機事件，可從任意分布進行抽樣，而不是嚴格按照物理概率抽樣，前提是需調整粒子權重來進行補償修正。因此，修正后的抽樣方法使得抽樣可按希望的方向發射粒子、粒子可進入感興趣的相空間（時間、能量），以此改變碰撞位置和類型。主要包括指數變換、隱俘獲、強迫碰撞、源偏倚和中子引發的光子產生偏倚。

4）部分確定論方法

該方法使用類似確定論方法來避開標準隨機游走過程。例如：下一事件估計、控制隨機數序列。MCNP 中該類方法包括點探測器、DXTRAN以及相關抽樣。

成功實現MCNP 的減方差方法通常是較為困難的，且藝術技巧性高于科學性［1］。每項降方差方法均有其自身的優點，但也有部分通用的原則。該通用原則就是要對物理問題有比較清晰的理解，依據物理本質來選取相應的方差降低方法。當選定了一項方法后，需指定適合的參數，這通常較選擇方法更困難。這些參數的首次猜測通常來自于相似問題的計算經驗，或同一問題的相似計算。指定好參數后，還需進行一系列短時間的運行以監測這些參數的有效性。

MCNP輸出文件包含了有助于了解抽樣的許多有用信息，需查看以下信息來檢測方差降低方法參數設置的有效性：改善了計數粒子的抽樣；各種降方差方法協作運行，不會嚴重相互影響；FOM 表穩定，該表可檢測抽樣不足；結果不會出現明顯錯誤。在確定降方差方法對計算問題能起到改善作用后，用戶還需進行多次短時間運行，以此逐步改善參數。這其中參數設置有效性的評判標準主要包括相對誤差R和品質因子FOM。

相對誤差R 是MCNP 計算結果可信度的評判標準，其定義為估計值的標準差與估計值的比值，即：

品質因子FOM 為MCNP 計算效率的評判標準，其定義為：

其中，T 為計算耗時。實際上，可從計數波動表（TFC）中監測相對誤差R 和品質因子FOM 隨粒子數的變化趨勢。

2.2 3種減方差方法

為獲取安全棒頂端的中子注量率情況，根據目前的建模方式，即整根棒為1個柵元，整個上方水池為1個柵元，選取3種減方差方法：點探測器、點探測器＋DXTRAN、點探測器＋強迫碰撞。

1）探測器貢獻［1］

由于粒子精確輸運到某點的概率是非常小的，因此使用偽粒子來定向到該點。每產生1個粒子（可以是源粒子，也可以是碰撞后的粒子），就要求有1個偽粒子在指定的空間位置點被計數。對小體積計數而言，很可能由于粒子到達該體積的概率極低而不能得到徑跡長度估計或穿過表面的通量估計，而點探測器仍可得到較好的通量估計，這就是其一大優點。正是由于每產生1個粒子就要在點探測器位置處進行計數操作，故將影響整體計算效率。

2）DXTRAN［1］

DXTRAN 即 確 定 性 輸 運（deterministic transport），它同時偏倚散射方向和源的方向。由于粒子幾乎散射不到一些小區域內，因此會造成小區域抽樣不充分。為改善這個情形，可在輸入文件中指定1 個DXTRAN 球，它需包含該小區域。一旦粒子在球外發生碰撞，或球外有源，DXTRAN 技巧就會生成1 個特殊的DXTRAN 粒子，且會確定性地散射到這個DXTRAN 球內。

3）強迫碰撞［1］

強迫碰撞方法將粒子分裂為碰撞部分和非碰撞部分。碰撞部分在當前柵元內發生強迫碰撞，非碰撞部分離開該柵元且不發生碰撞，但該部分是暫時存儲起來，等到該粒子的徑跡在柵元邊界上時才繼續追蹤。非碰撞部分繼續游走時的權為W0e-Σtd，其中W0為粒子初始權，Σt為宏觀總截面，d 為距離。即非碰撞部分的權等于當前粒子權乘以無碰撞地離開該柵元的概率。相應地，強迫碰撞部分的權為W0（1-e-Σtd），即總的權是守恒的。

3 計算結果分析

3.1 初步對比分析

為獲取安全棒最頂端的中子注量率情況，首先用較少的粒子數計算，檢驗3種方法的優劣。主要考察3種方法在相同計算條件下的相對誤差R 和品質因子FOM 隨粒子數nps的變化趨勢。計算方式為基于MPICH 的消息傳遞并行計算［4］，計算條件為：1）相同的臨界計算控制卡kcode 3000 1.0 10 40；2）相同的初始源，均來自同一個srctp源分布文件；3）相同的隨機數序列和隨機數種子；4）相同的柵元重要性，均設為1。

計算結果如圖3和4所示。圖中，F5代表點探測器計數；F5＋DXTRAN 代表點探測器加DXTRAN 減方差方法；F5＋FCL 代表點探測器加強迫碰撞減方差方法。

圖3 3種減方差方法的計算相對誤差Fig.3 Relative error of three variance reduction methods

圖4 3種減方差方法的FOMFig.4 FOM of three variance reduction methods

由圖3和4的計算結果可知，在相同的計算條件下，點探測器加強迫碰撞的相對誤差R較小，品質因子FOM 較高，其減方差效果較好。故采用該減方差組合進行增加粒子數計算，計算結果的相對誤差及品質因子如圖5所示。

圖5 增加粒子數運行結果Fig.5 Run result of increasing particle histories

由圖5可知，即便是在增加粒子數后，相對誤差R 也未明顯減小，呈現波動現象，同時，品質因子FOM 逐漸下降，即計算效率越來越低。且在粒子數增加到一定程度后，相對誤差發生階躍升高，計算獲得的結果數據發散。

由此可知，針對本實際問題，安全棒最頂端位置采用點探測器＋強迫碰撞的方式并不能達到預期的減方差效果，故需要采取其他方法，如賭分裂，將粒子引至感興趣的區域，但由于當前的幾何模型中安全棒是由1個柵元定義，且堆芯上部水層也是由1個柵元定義，直接指定柵元重要性并不能實現減方差效果，故考慮將其分層建模。

3.2 分層建模

通過重新建模，將堆芯上部水反射層分層建模，并通過重要性設置，希望進入頂層水的中子徑跡數增多［5］。這種方法的優點是，當劃分的區域和確定的重要性合適時，可使更多的粒子被引到感興趣區域的附近，從而可緩解常規隨機游動的粒子很少能到達計數區域附近的問題。為便于查看減方差效果，將安全棒頂端的小塊按樣品研究所需的大小（φ17.5 mm×30mm）單獨建模，稱該柵元為目標柵元。水反射層分層建模如圖6所示。

柵元重要性設置、中子在各水層和目標柵元中的徑跡及碰撞情況列于表1。第1 次計算時，將遠離堆芯的水層柵元重要性依次加倍，并將目標柵元與其所在水層的重要性均設置為最高的131 072。計算結果顯示各水層中的中子徑跡與碰撞情況都得到了提升，但目標柵元中的中子徑跡卻較少。為實現目標柵元中的中子徑跡與其同高度的水層相當，將目標柵元重要性調高為149 000／70×131 072＝2.79×108。從新的重要性分配下的第2次計算結果可知，目標柵元中的中子徑跡及碰撞數均達到較高值，起到了將中子引至目標柵元的效果。

再增加粒子數計算，以臨界控制條件kcode 3000 1.0 10 6010進行計算，對目標柵元采用徑跡長度體通量計數（F4），得到的計算結果相對誤差及品質因子如圖7所示。該體通量計數結果通過全部10個判斷標準，并從圖7可知：相對誤差隨運行粒子數的增加逐漸減小，且減小趨勢符合；品質因子FOM 雖整體較小，計算效率相對較低，但在收斂后基本保持穩定，也說明了該計算的穩定性。總體來講，依靠分層建模，進行柵元重要性分配，將粒子引向目標柵元，實現了減方差的目的。

圖6 水反射層分層建模Fig.6 Multiple layers modeling of water reflector

表1 柵元中的中子徑跡情況Table 1 Neutron track activity in cells

圖7 分層建模增加粒子數計算結果Fig.7 Run result of multiple layers modeling after increasing particle histories

4 結論

根據本研究堆幾何模型特點，首先，選取了點探測器、點探測器＋DXTRAN、點探測器＋強迫碰撞3種減方差方法進行對比計算，并對其中效果較好的點探測器＋強迫碰撞方法進行了進一步的增加粒子數計算，計算結果顯示：隨著粒子數的增加，計算發散，未起到減方差的作用。之后，考慮使用柵元重要性分配，將中子引向目標柵元。為實現此目的，將原有幾何模型重新分層建模，并分配適當的柵元重要性。計數結果表明：計數相對誤差在5%以內，品質因子保持穩定，實現了減方差的目的。

［1］ JUDITH F B.MCNP：A general Monte Carlo N-particle transport code version 4C，LA-13709-M［R］.USA：Los Alamos National Laboratory，2000.

［2］ 張嶄，李君利，武禎，等.內照射小器官劑量計算中減方差技巧的比較和應用［J］.清華大學學報：自然科學版，2007，47（S1）：1 051-1 056.ZHANG Zhan，LI Junli，WU Zhen，et al.Comparison and application of variance-reduction techniques used in internal radiation dose calculations for small organs［J］.Journal of Tsinghua University：Science and Technology，2007，47（S1）：1 051-1 056（in Chinese）.

［3］ 竇海峰，代君龍.SPRR-300 反應堆輻照孔道中子注量率的MCNP 程序計算［J］.核動力工程，2006，27（1）：51-54.DOU Haifeng，DAI Junlong.Calculation of neutron flux in SPRR-300reactor irradiation channel with MCNP code［J］.Nuclear Power Engineering，2006，27（1）：51-54（in Chinese）.

［4］ 楊萬奎，曾和榮，冷軍，等.300＃研究堆首爐中央孔道中子通量密度計算［J］.強激光與粒子束，2012，24（12）：3 000-3 005.YANG Wankui，ZENG Herong，LENG Jun，et al.Neutron flux calculation for central channel in first cycle of SPRR-300［J］.High Power Laser and Particle Beams，2012，24（12）：3 000-3 005（in Chinese）.

［5］ THOMAS E B.A sample problem for variance reduction in MCNP，LA-10363-MS［R］.USA：Los Alamos National Laboratory，1985.