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

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

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

1 堆芯布置及安全棒位置

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

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

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

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

2 減方差原理

2.1 一般原理

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

1）截?cái)喾椒?/p>

該類方法是最簡單直接的方法，主要包括幾何截?cái)?、能量截?cái)嗪蜁r間截?cái)唷?/p>

2）數(shù)量控制方法

該類方法通過粒子分裂和俄羅斯輪盤賭來控制抽樣數(shù)量。在重要的區(qū)域增加抽樣，但每個粒子的權(quán)重降低；在不重要的區(qū)域減少抽樣，但每個粒子的權(quán)重增大。該方法通過調(diào)整權(quán)重以保證結(jié)果無偏。主要包括幾何、能量、時間相關(guān)的賭分裂、權(quán)重截?cái)嘁约皺?quán)窗。

3）修正抽樣方法

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

4）部分確定論方法

該方法使用類似確定論方法來避開標(biāo)準(zhǔn)隨機(jī)游走過程。例如：下一事件估計(jì)、控制隨機(jī)數(shù)序列。MCNP 中該類方法包括點(diǎn)探測器、DXTRAN以及相關(guān)抽樣。

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

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

相對誤差R 是MCNP 計(jì)算結(jié)果可信度的評判標(biāo)準(zhǔn)，其定義為估計(jì)值的標(biāo)準(zhǔn)差與估計(jì)值的比值，即：

品質(zhì)因子FOM 為MCNP 計(jì)算效率的評判標(biāo)準(zhǔn)，其定義為：

其中，T 為計(jì)算耗時。實(shí)際上，可從計(jì)數(shù)波動表（TFC）中監(jiān)測相對誤差R 和品質(zhì)因子FOM 隨粒子數(shù)的變化趨勢。

2.2 3種減方差方法

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

1）探測器貢獻(xiàn)［1］

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

2）DXTRAN［1］

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

3）強(qiáng)迫碰撞［1］

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

3 計(jì)算結(jié)果分析

3.1 初步對比分析

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

計(jì)算結(jié)果如圖3和4所示。圖中，F(xiàn)5代表點(diǎn)探測器計(jì)數(shù)；F5＋DXTRAN 代表點(diǎn)探測器加DXTRAN 減方差方法；F5＋FCL 代表點(diǎn)探測器加強(qiáng)迫碰撞減方差方法。

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

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

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

圖5 增加粒子數(shù)運(yùn)行結(jié)果Fig.5 Run result of increasing particle histories

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

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

3.2 分層建模

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

柵元重要性設(shè)置、中子在各水層和目標(biāo)柵元中的徑跡及碰撞情況列于表1。第1 次計(jì)算時，將遠(yuǎn)離堆芯的水層?xùn)旁匾砸来渭颖?，并將目?biāo)柵元與其所在水層的重要性均設(shè)置為最高的131 072。計(jì)算結(jié)果顯示各水層中的中子徑跡與碰撞情況都得到了提升，但目標(biāo)柵元中的中子徑跡卻較少。為實(shí)現(xiàn)目標(biāo)柵元中的中子徑跡與其同高度的水層相當(dāng)，將目標(biāo)柵元重要性調(diào)高為149 000／70×131 072＝2.79×108。從新的重要性分配下的第2次計(jì)算結(jié)果可知，目標(biāo)柵元中的中子徑跡及碰撞數(shù)均達(dá)到較高值，起到了將中子引至目標(biāo)柵元的效果。

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

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

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

圖7 分層建模增加粒子數(shù)計(jì)算結(jié)果Fig.7 Run result of multiple layers modeling after increasing particle histories

4 結(jié)論

根據(jù)本研究堆幾何模型特點(diǎn)，首先，選取了點(diǎn)探測器、點(diǎn)探測器＋DXTRAN、點(diǎn)探測器＋強(qiáng)迫碰撞3種減方差方法進(jìn)行對比計(jì)算，并對其中效果較好的點(diǎn)探測器＋強(qiáng)迫碰撞方法進(jìn)行了進(jìn)一步的增加粒子數(shù)計(jì)算，計(jì)算結(jié)果顯示：隨著粒子數(shù)的增加，計(jì)算發(fā)散，未起到減方差的作用。之后，考慮使用柵元重要性分配，將中子引向目標(biāo)柵元。為實(shí)現(xiàn)此目的，將原有幾何模型重新分層建模，并分配適當(dāng)?shù)臇旁匾?。?jì)數(shù)結(jié)果表明：計(jì)數(shù)相對誤差在5%以內(nèi)，品質(zhì)因子保持穩(wěn)定，實(shí)現(xiàn)了減方差的目的。

［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］ 張嶄，李君利，武禎，等.內(nèi)照射小器官劑量計(jì)算中減方差技巧的比較和應(yīng)用［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，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 反應(yīng)堆輻照孔道中子注量率的MCNP 程序計(jì)算［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ì)算［J］.強(qiáng)激光與粒子束，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.