基于連續能量蒙特卡羅方法的均勻化群常數計算

李滿倉，王 侃，姚 棟

（1.清華大學工程物理系，北京100084；2.中國核動力研究設計院核反應堆系統設計技術重點實驗室，四川 成都610041）

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

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

以下首先介紹本文蒙特卡羅均勻化群常數的計算方法，之后在組件和堆芯兩個層次對方法進行驗證。

1 計算方法

1.1 一般群常數計算方法

本文一般群截面和群常數包括總截面、吸收截面、中子產生截面、裂變截面和中子裂變譜等，基本計算思路是使用徑跡長度方法估計反應率和通量，通過二者之商得到截面。

定義時間積分體積權重的群通量φg，

本文的所有計算都是靜態的，和時間無關。式（1）以及下面公式中的時間變量t，指的是粒子從產生到消失的時間。其中通量

改寫成距離的形式

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

其中，WTLiV（E）為第i個粒子在能量E體積V處的粒子權重和徑跡長度的乘積；Wi0是第i個粒子的初始權重；N為跟蹤的粒子總數。

類似的，對于反應率，有下式成立：

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

1.2 散射事件方法

群常數中的群間轉移截面和散射截面的高階勒讓德分量不能直接應用確定論的方法獲得。在連續能量點截面數據庫中沒有相應的微分散射截面，散射反應是通過總散射截面和散射角分布等量體現的；非彈性散射中的出射中子能量也有一定的分布。為此，論文提出散射事件方法計算散射矩陣和勒讓德分量。散射事件方法利用蒙特卡羅直接模擬的特點，跟蹤模擬散射事件，當散射發生時，統計能量和角度的變化，在此基礎上獲得散射矩陣和勒讓德系數。

對于散射矩陣，分別記錄散射事件前后的中子能量，判斷此能量所屬的分群能量區間并進行統計，獲得散射份額

群間轉移截面通過散射份額與散射截面的乘積獲得：

應用蒙特卡羅方法產生組件均勻化群常數時同時考慮下游堆芯應用擴散方法和輸運方法計算，因此有必要產生高階勒讓德分量，或稱Pn截面。

以一維問題為例，勒讓德分量的定義如下：

Σs，n和φn是使用勒讓德多項式展開散射截面和中子通量的系數，可以利用勒讓德多項式的正交性求得。其中Σs，0是群間轉移截面，φ0是中子注量率。n≥1時，應用蒙特卡羅方法統計權重函數φn（z，E）時，由于蒙特卡羅方法的漲落特性，會造成很大誤差。假設φn與φ0成比例關系，即

應用微分散射截面的定義

（10）式可以寫成

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

1.3 擴散系數

擴散系數是輸運方程到擴散方程近似時引入的參數。擴散系數定義為：

其中Σtr，g為輸運截面。確定論方法使用多群核數據庫中包含有各核素的輸運截面，而在連續能量蒙特卡羅數據庫中沒有對應的量。作為初步實現蒙特卡羅方法組件均勻化的方法，本文輸運截面通過下式獲得：

其中 Σt，g為總截面，Σs，1，g→g是P1系數的自散射截面，應用1.2節所述計算方法獲得。

1.4 泄漏修正

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

多群形式的B1方程

其中

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

臨界譜是考慮了臨界效應的能譜，應用該能譜作為群常數歸并時的權重函數：

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

2 數值驗證

2.1 驗證方法

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

MCMC具有產生任意數目群參數的能力。表1給出了2群、4群和7群的能群邊界，本文將考察3種能群結構數對物理計算結果的影響。

表1 2群、4群和7群的能群結構Table1 Group boundaries for the 2－，4－and 7－group structures

組件計算時如果應用全反射邊界條件，凈流為零，中子輸運方程和擴散方程中的泄漏項為零。以多群擴散方程為例：）

對于全反射邊界條件，方程（20）寫作

即

根據（22）式可以根據均勻化群常數計算kinf的理論值。這可以作為另一個驗證組件群常數有效性的手段。

2.2 AFA3G組件驗證

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

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

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

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

表3給出了CASMO－3和MCMC計算的群常數解析解kinf與參考值的比較。采用MCMC群常數計算的解析解kinf與參考值符合較好，隨著能群數目的增加，結果越來越好；總體來說，CASMO群常數的結果要差一些，而且沒有特征值隨能群數目增加而改善的趨勢。這是因為確定論組件計算時采用多群數據庫，MCMC計算應用連續能量，均勻化歸并少群參數時保留了更多連續能量的信息。

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

2.3 組件層次幾何適應性研究

圖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少群常數的kinf解析解和參考解的比較在表4中給出。結果表明：對不同幾何構型和材料組成的組件，MCMC群常數反算kinf與MCNP參考值符合均較好，隨能群數目增加而更加吻合。這說明MCMC計算的組件群常數（總截面、吸收截面、中子產生截面、散射矩陣和裂變譜）是準確有效的。相比目前實際工作中每種燃料組件的計算都需要專門的組件程序而言，MCMC有良好的幾何適應性。

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

2.4 簡化壓水堆堆芯驗證

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

圖5 簡化壓水堆模型堆芯構成Fig.5 Simplified PWR core configuration

表5是堆芯keff的計算結果。表中MCMC群常數和CASMO群常數都是經過泄漏修正的，MCMC群常數沒有應用任何形式的等效均勻化，CASMO群常數使用G因子方法調整相對吸收率守恒[20］。MCMG計算中比較加入了P1階勒讓德分量的效果。2、4和7表示使用的少群群常數分別是2群、4群和7群。可以看到，應用MCMC群常數的輸運計算和擴散計算keff基本隨能群數目增加與MCNP參考值的相對偏差變小。應用MCMC群常數的擴散計算精度優于確定論CASMO程序。加入P1高階分量的MCMG計算堆芯keff比P0的計算結果好很多，因為P1截面的加入會改善角度信息。

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

表5 簡化壓水堆堆芯keff比較Table5 keffcomparison for simplified PWR core

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

3 結束語

本文旨在將蒙特卡羅方法對復雜幾何和物理的適應性以及連續能量的優勢應用于傳統的兩步法反應堆物理計算流程中，根據蒙特卡羅方法的特點總結和提出了蒙特卡羅方法產生組件均勻化群常數的方法。通過組件和堆芯兩個尺度的數值驗證，表明本文計算的蒙特卡羅組件群常數具有較好的性能，在幾何上有很強的適應性，而且保留了更多連續能量信息。基于這些優勢，蒙特卡羅方法均勻化和確定論組件堆芯計算的反應堆物理分析流程，是下一代堆芯設計計算中極有潛力的候選方法。

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

[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中子輸運蒙特卡羅程序 MCMG－I使用說明[R］.北京：北京應用物理與計算數學研究所，2011.

[18]張瑞茵，趙翊民.CITATION核反應堆堆芯分析程序[R］.清華大學核能與新能源研究院，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.