李滿倉(cāng)，王 侃，姚 棟

（1.清華大學(xué)工程物理系，北京100084；2.中國(guó)核動(dòng)力研究設(shè)計(jì)院核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都610041）

在反應(yīng)堆物理計(jì)算中應(yīng)用的確定論方法和蒙特卡羅方法方法各有優(yōu)缺點(diǎn)：確定論方法一般計(jì)算速度較快，蒙特卡羅方法通常精度上具有優(yōu)勢(shì)。傳統(tǒng)的反應(yīng)堆物理分析以確定論為主，蒙特卡羅方法通常只是作為補(bǔ)充手段。隨著計(jì)算機(jī)技術(shù)和高性能算法的不斷發(fā)展，完成物理計(jì)算所需時(shí)間越來(lái)越短，人們希望在更高的平臺(tái)上追求速度與精度的平衡；另一方面，新型復(fù)雜堆芯概念的提出，對(duì)程序所能模擬的幾何和物理復(fù)雜度也形成挑戰(zhàn)。這種背景下，聯(lián)合蒙特卡羅方法和確定論方法優(yōu)勢(shì)的研究是近幾年的一個(gè)熱點(diǎn)方向[1－3］。其中一個(gè)思路是應(yīng)用蒙特卡羅方法計(jì)算組件為確定論堆芯程序提供均勻化群常數(shù)。這種方案沿用傳統(tǒng)的“兩步法”，提高對(duì)復(fù)雜幾何和物理的適應(yīng)性，避免復(fù)雜的共振計(jì)算，同時(shí)保證堆芯計(jì)算在可接受的時(shí)間內(nèi)完成。使用同一個(gè)蒙特卡羅組件均勻化程序和同一套數(shù)據(jù)即可模擬計(jì)算任意幾何構(gòu)型的組件和任意能譜的堆芯，普適性很好。

近幾年國(guó)際上一直在探索應(yīng)用蒙特卡羅方法進(jìn)行組件均勻化[4－11］，應(yīng)用涉及沸水堆、高通量堆、乏燃料組件和球床堆等。原則上所有的蒙特卡羅程序都可以通過(guò)計(jì)數(shù)能量和體積區(qū)間的反應(yīng)率和通量獲得均勻化截面。但是廣泛使用的蒙特卡羅程序如 MCNP[12］、KENO[13］、TRIPOLI[14］、MVP[15］等不具備這樣的功能。應(yīng)用連續(xù)能量蒙特卡羅方法進(jìn)行組件尺度的均勻化計(jì)算主要涉及兩個(gè)大的問(wèn)題：一是群常數(shù)的基本產(chǎn)生方法，蒙特卡羅方法相對(duì)于確定論有一定的特殊性，針對(duì)確定論發(fā)展的均勻化理論并不能直接應(yīng)用在蒙特卡羅方法中；二是群常數(shù)等效均勻化，均勻化過(guò)程必須做等效處理，蒙特卡羅方法均勻化也不例外。本文研究和探索解決第一個(gè)問(wèn)題，即蒙特卡羅方法計(jì)算產(chǎn)生組件均勻化群常數(shù)，特別是散射矩陣和高階勒讓德系數(shù)。擴(kuò)散系數(shù)在蒙特卡羅方法中沒(méi)有對(duì)應(yīng)的物理量，通過(guò)研究采用P1截面計(jì)算。對(duì)于全反射單組件計(jì)算模型，本文應(yīng)用BN理論對(duì)均勻化群常數(shù)進(jìn)行泄漏修正。

以下首先介紹本文蒙特卡羅均勻化群常數(shù)的計(jì)算方法，之后在組件和堆芯兩個(gè)層次對(duì)方法進(jìn)行驗(yàn)證。

1 計(jì)算方法

1.1 一般群常數(shù)計(jì)算方法

本文一般群截面和群常數(shù)包括總截面、吸收截面、中子產(chǎn)生截面、裂變截面和中子裂變譜等，基本計(jì)算思路是使用徑跡長(zhǎng)度方法估計(jì)反應(yīng)率和通量，通過(guò)二者之商得到截面。

定義時(shí)間積分體積權(quán)重的群通量φg，

本文的所有計(jì)算都是靜態(tài)的，和時(shí)間無(wú)關(guān)。式（1）以及下面公式中的時(shí)間變量t，指的是粒子從產(chǎn)生到消失的時(shí)間。其中通量

改寫(xiě)成距離的形式

其中s是粒子從產(chǎn)生到消失的距離。在蒙特卡羅方法中，粒子密度等價(jià)于每單位體積的粒子權(quán)重W的總和，而距離等價(jià)于徑跡長(zhǎng)度TL。因此（3）式可以寫(xiě)成

其中，WTLiV（E）為第i個(gè)粒子在能量E體積V處的粒子權(quán)重和徑跡長(zhǎng)度的乘積；Wi0是第i個(gè)粒子的初始權(quán)重；N為跟蹤的粒子總數(shù)。

類似的，對(duì)于反應(yīng)率，有下式成立：

根據(jù)（5）式和（6）式即可得到均勻化群截面

1.2 散射事件方法

群常數(shù)中的群間轉(zhuǎn)移截面和散射截面的高階勒讓德分量不能直接應(yīng)用確定論的方法獲得。在連續(xù)能量點(diǎn)截面數(shù)據(jù)庫(kù)中沒(méi)有相應(yīng)的微分散射截面，散射反應(yīng)是通過(guò)總散射截面和散射角分布等量體現(xiàn)的；非彈性散射中的出射中子能量也有一定的分布。為此，論文提出散射事件方法計(jì)算散射矩陣和勒讓德分量。散射事件方法利用蒙特卡羅直接模擬的特點(diǎn)，跟蹤模擬散射事件，當(dāng)散射發(fā)生時(shí)，統(tǒng)計(jì)能量和角度的變化，在此基礎(chǔ)上獲得散射矩陣和勒讓德系數(shù)。

對(duì)于散射矩陣，分別記錄散射事件前后的中子能量，判斷此能量所屬的分群能量區(qū)間并進(jìn)行統(tǒng)計(jì)，獲得散射份額

群間轉(zhuǎn)移截面通過(guò)散射份額與散射截面的乘積獲得：

應(yīng)用蒙特卡羅方法產(chǎn)生組件均勻化群常數(shù)時(shí)同時(shí)考慮下游堆芯應(yīng)用擴(kuò)散方法和輸運(yùn)方法計(jì)算，因此有必要產(chǎn)生高階勒讓德分量，或稱Pn截面。

以一維問(wèn)題為例，勒讓德分量的定義如下：

Σs，n和φn是使用勒讓德多項(xiàng)式展開(kāi)散射截面和中子通量的系數(shù)，可以利用勒讓德多項(xiàng)式的正交性求得。其中Σs，0是群間轉(zhuǎn)移截面，φ0是中子注量率。n≥1時(shí)，應(yīng)用蒙特卡羅方法統(tǒng)計(jì)權(quán)重函數(shù)φn（z，E）時(shí)，由于蒙特卡羅方法的漲落特性，會(huì)造成很大誤差。假設(shè)φn與φ0成比例關(guān)系，即

應(yīng)用微分散射截面的定義

（10）式可以寫(xiě)成

（13）式就是本文Pn的計(jì)算方法：首先在－1到1之間分割一定數(shù)目的角度區(qū)間，根據(jù)散射事件前后角度的改變獲得散射角余弦，在相應(yīng)的角度區(qū)間上增加計(jì)數(shù)；粒子跟蹤完畢，可以獲得各個(gè)角度區(qū)間的份額，然后利用勒讓德關(guān)系式的定義，計(jì)算數(shù)值積分。

1.3 擴(kuò)散系數(shù)

擴(kuò)散系數(shù)是輸運(yùn)方程到擴(kuò)散方程近似時(shí)引入的參數(shù)。擴(kuò)散系數(shù)定義為：

其中Σtr，g為輸運(yùn)截面。確定論方法使用多群核數(shù)據(jù)庫(kù)中包含有各核素的輸運(yùn)截面，而在連續(xù)能量蒙特卡羅數(shù)據(jù)庫(kù)中沒(méi)有對(duì)應(yīng)的量。作為初步實(shí)現(xiàn)蒙特卡羅方法組件均勻化的方法，本文輸運(yùn)截面通過(guò)下式獲得：

其中 Σt，g為總截面，Σs，1，g→g是P1系數(shù)的自散射截面，應(yīng)用1.2節(jié)所述計(jì)算方法獲得。

1.4 泄漏修正

組件均勻化計(jì)算一般采用全反射單組件局部均勻化模型，這種模型沒(méi)有考慮組件在堆芯中的真實(shí)狀態(tài)：組件的邊界存在中子凈流，組件通常處于臨界狀態(tài)。蒙特卡羅方法均勻化也要還原兩步法計(jì)算流程中組件在堆芯的臨界狀態(tài)，本文采用確定論方法中成熟的BN理論修正泄漏效應(yīng)，基本思路是：使用B1方程中的曲率B2表征堆芯的有限性，使用近似漸近能譜考慮中子泄漏。

多群形式的B1方程

其中

（16）式和（17）式中B2為臨界曲率，其他參量與一般多群擴(kuò)散方程中的意義相同。通過(guò)求解（16）式獲得臨界漸近能譜φLn，使用臨界漸進(jìn)譜修正無(wú)限介質(zhì)譜，得到臨界譜φCn，i：

臨界譜是考慮了臨界效應(yīng)的能譜，應(yīng)用該能譜作為群常數(shù)歸并時(shí)的權(quán)重函數(shù)：

這樣在全反射單組件模型下得到的均勻化群常數(shù)就可以歸并出考慮了臨界效應(yīng)的少群常數(shù)。該少群常數(shù)用于堆芯計(jì)算即可提高其精度。

2 數(shù)值驗(yàn)證

2.1 驗(yàn)證方法

數(shù)值驗(yàn)證部分蒙特卡羅方法均勻化群常數(shù)采用MCMC程序計(jì)算。MCMC程序是本文作者編制的基于連續(xù)能量蒙特卡羅方法組件均勻化程序。為展示本文研究的方法和蒙特卡羅均勻化的特性，以下計(jì)算中MCMC群常數(shù)沒(méi)有應(yīng)用等效均勻化。MCMC計(jì)算獲得的均勻化群常數(shù)將與CASMO－3[16］進(jìn)行對(duì)比。應(yīng)用多群蒙特卡羅程序MCMG[17］和細(xì)網(wǎng)有限差分程序CITATION[18］計(jì)算堆芯，獲得keff和通量分布，與參考解MCNP[19］的計(jì)算結(jié)果對(duì)比，驗(yàn)證方法的可靠性。

MCMC具有產(chǎn)生任意數(shù)目群參數(shù)的能力。表1給出了2群、4群和7群的能群邊界，本文將考察3種能群結(jié)構(gòu)數(shù)對(duì)物理計(jì)算結(jié)果的影響。

表1 2群、4群和7群的能群結(jié)構(gòu)Table1 Group boundaries for the 2－，4－and 7－group structures

組件計(jì)算時(shí)如果應(yīng)用全反射邊界條件，凈流為零，中子輸運(yùn)方程和擴(kuò)散方程中的泄漏項(xiàng)為零。以多群擴(kuò)散方程為例：）

對(duì)于全反射邊界條件，方程（20）寫(xiě)作

即

根據(jù)（22）式可以根據(jù)均勻化群常數(shù)計(jì)算kinf的理論值。這可以作為另一個(gè)驗(yàn)證組件群常數(shù)有效性的手段。

2.2 AFA3G組件驗(yàn)證

圖1所示為17×17形式的AFA3G組件，富集度為2%。分別使用CASMO－3和MCMC計(jì)算產(chǎn)生控制棒全提（ARO）和控制棒全插（ARI）兩種工況的均勻化群常數(shù)。

圖1 AFA3G組件模型Fig.1 AFA3Gassembly model

表2比較了MCMC與CASMO計(jì)算的2群群常數(shù)。從表中可以看到：無(wú)論是ARO工況還是ARI工況，除了跳群散射截面，MCMC與CASMO－3計(jì)算的群常數(shù)相對(duì)偏差都在5%以內(nèi)，對(duì)于截面而言，這是很小的；不同群的群間轉(zhuǎn)移截面差別較大，考慮到熱群到快群向上散射截面本身數(shù)值就很小，在很多2群堆芯計(jì)算中甚至設(shè)為0，其影響并不顯著。

表2 AFA3G組件CASMO－3和MCMC計(jì)算的2群群常數(shù)比較Table2 Comparison between two groups of constants in MCMC and CASMO－3calculations for AFA3Gassembly

表3給出了CASMO－3和MCMC計(jì)算的群常數(shù)解析解kinf與參考值的比較。采用MCMC群常數(shù)計(jì)算的解析解kinf與參考值符合較好，隨著能群數(shù)目的增加，結(jié)果越來(lái)越好；總體來(lái)說(shuō)，CASMO群常數(shù)的結(jié)果要差一些，而且沒(méi)有特征值隨能群數(shù)目增加而改善的趨勢(shì)。這是因?yàn)榇_定論組件計(jì)算時(shí)采用多群數(shù)據(jù)庫(kù)，MCMC計(jì)算應(yīng)用連續(xù)能量，均勻化歸并少群參數(shù)時(shí)保留了更多連續(xù)能量的信息。

表3 AFA3G組件kinf解析解與參考值比較Table3 Theoretical solutions ofkinffor AFA3Gassembly

2.3 組件層次幾何適應(yīng)性研究

圖2是套筒狀燃料元件模型，燃料和包殼呈套筒狀相間排列；圖3為CANFLEX重水堆新型燃料棒束組件模型，含43根燃料棒；圖4所示的雙排六邊形燃料組件是超臨界水堆的候選組件。

圖2 套筒狀燃料元件模型Fig.2 Cylindrical assembly model

圖3 CANFLEX組件模型Fig.3 CANFLEX assembly model

圖4 SCWR組件模型Fig.4 SCWR assembly model

3種組件MCMC少群常數(shù)的kinf解析解和參考解的比較在表4中給出。結(jié)果表明：對(duì)不同幾何構(gòu)型和材料組成的組件，MCMC群常數(shù)反算kinf與MCNP參考值符合均較好，隨能群數(shù)目增加而更加吻合。這說(shuō)明MCMC計(jì)算的組件群常數(shù)（總截面、吸收截面、中子產(chǎn)生截面、散射矩陣和裂變譜）是準(zhǔn)確有效的。相比目前實(shí)際工作中每種燃料組件的計(jì)算都需要專門(mén)的組件程序而言，MCMC有良好的幾何適應(yīng)性。

表4 套筒狀燃料元件、CANFLEX組件和SCWR組件理論kinfTable4 Theoretical solutions of kinf for cylindrical assembly，CANFLEX and SCWR assembly

2.4 簡(jiǎn)化壓水堆堆芯驗(yàn)證

堆芯尺度驗(yàn)證在圖5所示的簡(jiǎn)化壓水堆模型上進(jìn)行。使用17×17形式的AFA3G組件，富集度為1.2%和2.0%的組件分別有24個(gè)和21個(gè)，兩種組件按棋盤(pán)布置。分別計(jì)算ARO和ARI兩種工況，ARI工況時(shí)2.0%組件中插入Ag－In－Cd控制棒。反射層為20cm厚的輕水。利用MCMC產(chǎn)生的群常數(shù)進(jìn)行多群蒙特卡羅MCMG和細(xì)網(wǎng)差分CITATION計(jì)算，比較堆芯keff和通量分布。反射層計(jì)算采用超組件模型。

圖5 簡(jiǎn)化壓水堆模型堆芯構(gòu)成Fig.5 Simplified PWR core configuration

表5是堆芯keff的計(jì)算結(jié)果。表中MCMC群常數(shù)和CASMO群常數(shù)都是經(jīng)過(guò)泄漏修正的，MCMC群常數(shù)沒(méi)有應(yīng)用任何形式的等效均勻化，CASMO群常數(shù)使用G因子方法調(diào)整相對(duì)吸收率守恒[20］。MCMG計(jì)算中比較加入了P1階勒讓德分量的效果。2、4和7表示使用的少群群常數(shù)分別是2群、4群和7群?？梢钥吹剑瑧?yīng)用MCMC群常數(shù)的輸運(yùn)計(jì)算和擴(kuò)散計(jì)算keff基本隨能群數(shù)目增加與MCNP參考值的相對(duì)偏差變小。應(yīng)用MCMC群常數(shù)的擴(kuò)散計(jì)算精度優(yōu)于確定論CASMO程序。加入P1高階分量的MCMG計(jì)算堆芯keff比P0的計(jì)算結(jié)果好很多，因?yàn)镻1截面的加入會(huì)改善角度信息。

表6是通量分布的偏差比較，包括最大偏差MAX、平均偏差A(yù)VG和均方根偏差RMS其中en為 MCMG或CITATION與MCNP參考值的相對(duì)誤差）。從結(jié)果看，CITATION與參考值偏差較小，而MCMG在少能群數(shù)目時(shí)偏差較大，需要更高階Pn截面完成堆芯輸運(yùn)計(jì)算。由于蒙特卡羅方法統(tǒng)計(jì)漲落的原因，MCNP參考值本身在堆芯對(duì)稱位置的分布也有約1%的偏差，所以總體而言，通量分布符合的很好。

表5 簡(jiǎn)化壓水堆堆芯keff比較Table5 keffcomparison for simplified PWR core

表6 簡(jiǎn)化壓水堆堆芯歸一化通量分布偏差比較Table6 Normalized flux distribution error comparison for simplified PWR core

3 結(jié)束語(yǔ)

本文旨在將蒙特卡羅方法對(duì)復(fù)雜幾何和物理的適應(yīng)性以及連續(xù)能量的優(yōu)勢(shì)應(yīng)用于傳統(tǒng)的兩步法反應(yīng)堆物理計(jì)算流程中，根據(jù)蒙特卡羅方法的特點(diǎn)總結(jié)和提出了蒙特卡羅方法產(chǎn)生組件均勻化群常數(shù)的方法。通過(guò)組件和堆芯兩個(gè)尺度的數(shù)值驗(yàn)證，表明本文計(jì)算的蒙特卡羅組件群常數(shù)具有較好的性能，在幾何上有很強(qiáng)的適應(yīng)性，而且保留了更多連續(xù)能量信息?；谶@些優(yōu)勢(shì)，蒙特卡羅方法均勻化和確定論組件堆芯計(jì)算的反應(yīng)堆物理分析流程，是下一代堆芯設(shè)計(jì)計(jì)算中極有潛力的候選方法。

連續(xù)能量蒙特卡羅方法均勻化產(chǎn)生群常數(shù)后，也要應(yīng)用等效均勻化理論保證均勻化前后的反應(yīng)率和界面流守恒，以最大限度的保留精度，體現(xiàn)蒙特卡羅方法的適用性。另外，均勻化群常數(shù)隨燃耗變化，為完成堆芯的燃耗計(jì)算，蒙特卡羅均勻化也要具備燃耗功能。這些是完善蒙特卡羅方法均勻化的研究方向。

[1]Brown F B，Martin W R，Mosteller R D.Monte Carlo－Advances and Challenges[C］／／Workshop at PHYSOR－2008，14－19September，Interlaken，Switzerland.Los Alamos National Laboratory，2008：159.

[2]Martin B.Advances in Monte Carlo Methods for Global Reactor Analysis [C］／／Joint International Topical Meeting on Mathematics ＆Computation and Supercomputing in Nuclear Applications（M＆C＋SNA 2007），April 15－19Monterey，California，USA.2007.

[3]Yesilyurt G，Martin W R，Lee J C.A Coupled Monte Carlo／Collision Probability Method for VHTR Analysis[J］.Transactions of the American Nuclear Society，2008，99：753－754.

[4]Ilas G，Rahnema F.A Monte Carlo Based Nodal Diffusion Model for Criticality Analysis of Spent Fuel Storage Lattices[J］.Annals of Nuclear Energy，2003，30（10）：1089－1108.

[5]Van der Marck S C，Kuijper J C，Oppe J.Homogenized Group Cross Section by Monte Carlo[C］／／PHYSOR－2006American Nuclear Society's Topical Meeting on Reactor Physics， Sept.10－14， Vancouver， BC，Canada.2006.

[6]Tohjoh M，Watanabe M，Yamamoto A.Application of Continuous－energy Monte Carlo Code as a Cross－section Generator of BWR Core Calculations[J］.Annals of Nuclear Energy，2005，32（8）：857－875.

[7]Redmond II E L.Multigroup Cross Section Generation Via Monte Carlo Methods[D］.Massachusetts Institute of Technology Dept.of Nuclear Engineering，1997.

[8]Rahnema F，Ilas G，Hudson N H，et al.On the Fewgroup Cross－section Generation Methodology for PBR Analysis[J］.Annals of Nuclear Energy，2006，33（11－12）：1058－1070.

[9]Blomquist R N.VIM Monte Carlo Neutron／Photon Transport Code[EB／OL］.[2010］.www.vim.anl.gov.

[10]Maiorov L V.MCU Code：MCU Project Home Page[EB／OL］.[2010］.www.mcu.vver.kiae.ru／eareal.

[11]Lepp?nen J.PSG2／Serpent－A Monte Carlo Reactor Physics Burnup Calculation Code[EB／OL］.[2010］.www.montecarlo.vtt.fi.

[12]Team X－M C.MCNP— A General Monte Carlo NParticle Transport Code， Version 5， Volume I：Overview and Theory[R］.LA－UR－03－1987，April 24，2003.

[13]Hollenbach D F，Petrie L M，Landers N F.KENOVI：A General Quadratic Version of the KENO Program[R］.ORNL／NUREG／CSD－2／V2／R6，Oak Ridge National Laboratory，1998.

[14]Diop C M，Petit O，Dumonteil E，et al.TRIPOLI－4：a 3DContinuous－Energy Monte Carlo Transport Code[C］／／PHYTRA1：First International Conference on Physics and Technology of Reactors and Applications，March 14－16，Marrakech，Morocco.2007.

[15]Matsuda S.MVP／GMVPⅡ：General Purpose Monte Carlo Codes for Neutron and Photon Transport Calculations based on Continuous Energy and Multigroup Methods[R］.JAERI 1348，Japan Atomic Energy Research Institute，2005.

[16]Studsvik.CASMO－3User's Manual[R］.Studsvik／SOA－94／9，Studsvik of America，1994.

[17]鄧力 .三維多群P5中子輸運(yùn)蒙特卡羅程序 MCMG－I使用說(shuō)明[R］.北京：北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，2011.

[18]張瑞茵，趙翊民.CITATION核反應(yīng)堆堆芯分析程序[R］.清華大學(xué)核能與新能源研究院，1986.

[19]Team X－M C.MCNP—A General Monte Carlo NParticle Transport Code，Version 5，Volume Ⅱ：User's Guide[R］.LA－UR－03－1987，April 24，2003.

[20]Studsvik.CASMO－3Methodology [R］.Studsvik／NFA－88／49，Studsvik，1988.