強(qiáng)勝龍，秦 冬，柴曉明，姚 棟

(中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都 610041)

用蒙特卡羅方法進(jìn)行輸運計算具有非常高的精度，但在耦合燃耗計算功能后，其計算的燃耗精度差強(qiáng)人意，特別是隨著燃耗的加深，誤差逐漸增大[1]。燃耗計算誤差影響核設(shè)計的精度，進(jìn)而影響整個堆芯的設(shè)計裕量。因此，在發(fā)展計算程序時，有必要對基于蒙特卡羅方法的燃耗(簡稱蒙卡燃耗)計算的誤差進(jìn)行深入分析。

蒙卡燃耗計算自20世紀(jì)80年代發(fā)展至今，誤差研究逐漸深入。Takeda等[2]提出燃耗矩陣法，指出蒙卡燃耗誤差的3個主要來源(截面、核密度、統(tǒng)計)并逐項分析，但該方法的數(shù)學(xué)推導(dǎo)較復(fù)雜，難以在程序中實現(xiàn)。Tohjoh等[3]提出了誤差的傳遞效應(yīng)，并給出估計誤差傳遞的公式，但未提出如何修正誤差。Dumonteil等[4]研究了蒙卡燃耗誤差中統(tǒng)計誤差的傳遞效應(yīng)，并針對蒙卡燃耗計算過程中隨機(jī)過程導(dǎo)致的偏差進(jìn)行了分析，選擇無偏差最小方差(UMV)算子用于求解波爾茲曼/貝曼方程。然而，Dumonteil等僅采用簡單的數(shù)學(xué)抽樣模型對UMV修正進(jìn)行模擬研究，并未在實際的蒙卡燃耗程序中進(jìn)行研究或應(yīng)用。

本工作參考Dumonteil等的工作，對蒙卡燃耗所獨有的誤差進(jìn)行研究，同時還針對具體的蒙卡燃耗計算程序提出相應(yīng)的修正措施。

1 蒙卡燃耗誤差修正研究

1.1 蒙卡燃耗誤差分類

蒙卡燃耗誤差的來源主要有3個：統(tǒng)計誤差、截面誤差和核密度誤差。這3類誤差相互影響，并隨燃耗計算而傳遞。從誤差的形成機(jī)制來看，蒙卡燃耗計算程序的誤差主要有3個層次：1) 方法層次，蒙特卡羅方法導(dǎo)致的統(tǒng)計誤差及其傳遞偏差；2) 數(shù)據(jù)庫層次，中子連續(xù)能譜點截面庫中裂變產(chǎn)額數(shù)據(jù)、單群截面數(shù)據(jù)等的誤差；3) 程序?qū)哟?，燃耗鏈方程求解及輸運燃耗方程求解的誤差。

其中，數(shù)據(jù)庫和方程求解的誤差不是蒙卡燃耗所獨有，確定論方法同樣存在該類問題。目前針對這兩類誤差的研究很多[5-7]。而統(tǒng)計誤差則是蒙特卡羅方法所獨有，特別是其傳遞效果。

1.2 UMV修正

若消除蒙卡燃耗計算過程中正態(tài)均值(通量φ)的非線性函數(shù)(核素含量NB(t))的偏差，可進(jìn)行UMV修正，重點是對UMV修正系數(shù)進(jìn)行計算。

本工作以MOI程序為基礎(chǔ)進(jìn)行研究。MOI是基于連續(xù)能量蒙特卡羅方法開發(fā)的堆芯調(diào)棒臨界燃耗計算程序，在UNIX或LINUX平臺下并行執(zhí)行，可實現(xiàn)調(diào)棒臨界燃耗計算。MOI計算的基本流程示于圖1。

圖1 MOI計算的基本流程

UMV修正系數(shù)包含σ、Ι、σ、Z、t，分別為轉(zhuǎn)移截面(包括中子吸收和散射)矩陣、中子的模擬代數(shù)、中子通量的標(biāo)準(zhǔn)方差、滿足(0,1)正態(tài)分布的隨機(jī)向量和燃耗步長。具體計算步驟如下：

1) 從MOI的燃耗計算輸入文件中獲得時間步長t。

2) 從MOI輸運計算輸入中讀取中子模擬代數(shù)I，通過蒙特卡羅模擬得到向量Z。

4) 在MOI燃耗計算輸出中讀取矩陣σ，求出UMV修正系數(shù)矩陣cos(σξt)。其中：

Z由正態(tài)分布中抽樣得到，正態(tài)分布示于圖2。

圖2 正態(tài)分布示意圖

5) MOI燃耗計算輸出的核素含量乘以修正系數(shù)(矩陣)，即可得到修正后的核素含量NNB(t)：

NNB(t)=NB(t)cos(σξt)

UMV修正流程示于圖3。

以一體式可燃毒物(IFBA)單柵元例題為例(圖4，UO2芯塊富集度為3.2%)，采用UMV修正后238U的相對修正量((修正后質(zhì)量-修正前質(zhì)量)/修正前質(zhì)量)隨燃耗的變化示于圖5，達(dá)到了Tohjoh等[3]觀察到的10-5量級。這說明UMV修正可改善核素含量。

圖3 UMV修正流程

圖4 IFBA柵元示意圖

圖5 238U修正量隨燃耗的變化

但UMV修正對kinf的影響甚微。仍以IFBA單柵元例題為例，與西屋公司開發(fā)的Paragon程序進(jìn)行對比驗證，同時還給出了阿海琺公司APOLLO2-F程序的計算結(jié)果(圖6)。Paragon與APOLLO2-F均廣泛應(yīng)用于商用核電廠的設(shè)計和后續(xù)技術(shù)服務(wù)，大量的運行和實驗數(shù)據(jù)反饋保證了程序的準(zhǔn)確性，因此選擇這兩個程序作為參考。計算結(jié)果顯示，采用UMV修正后的kinf仍存在較大誤差。

1.3 密度修正

蒙卡燃耗計算程序除統(tǒng)計誤差及其傳遞效應(yīng)外，其燃耗核素的選擇也會影響程序的精度。不同程序采取的策略不同：如MCODE[8]采取固定的核素選??；MCBurn[9]采用核素篩選三原則(優(yōu)先考慮吸收率，其次質(zhì)量分?jǐn)?shù)，以及重要核素等)；RMC[10]采用按宏觀吸收率份額篩選。

圖6 kinf隨燃耗的變化

無論采用何種策略，燃耗核素越多，則燃耗計算精度越高，但相應(yīng)的計算效率越低。以MCODE為例，在計算某單柵元例題時，為了高精度，選擇了39種錒系核素和100種裂變產(chǎn)物，保證了99.95%的質(zhì)量和中子吸收份額[8]。然而，該燃耗計算耗時4 d(DEC alpha工作站，每次MCNP計算耗時1 h)，在處理較大規(guī)模的問題時，該篩選策略會導(dǎo)致計算時間無法接受。

圖7 燃料芯塊相對質(zhì)量隨燃耗的變化

MCBurn等程序優(yōu)先考慮以吸收率份額篩選，在一定程度上減少了不必要的核素計算。然而，該篩選策略可能導(dǎo)致不參與輸運計算的裂變產(chǎn)物逐漸增多，因此篩選后的燃料質(zhì)量隨燃耗的增加而降低(圖7，按99.99%吸收率份額進(jìn)行篩選)。蒙特卡羅輸運計算通常輸入核素的相對含量，一旦總質(zhì)量發(fā)生變化，必須改變相應(yīng)燃料區(qū)的密度，才能真實反應(yīng)核素的質(zhì)量。

MOI采用按吸收率份額篩選的原則(同時也考慮重要核素如241Am、243Am等，該類核素隨燃耗逐漸積累，在深燃耗下作用較大)，因此隨燃耗的變化，MOI按照各燃耗區(qū)核素質(zhì)量的不斷變化修正各燃耗區(qū)的密度(簡稱密度修正)，從而保證各燃耗步下蒙特卡羅輸運計算的準(zhǔn)確性。

密度修正的公式為：

通過修正各燃耗區(qū)的密度，可避免密度誤差導(dǎo)致核素質(zhì)量的輸入誤差，保證蒙特卡羅輸運計算的準(zhǔn)確性。

1.4 誤差修正效果

通過對蒙卡燃耗統(tǒng)計誤差傳遞效應(yīng)的研究，提出了改善核密度的UMV修正；同時針對核素篩選導(dǎo)致的質(zhì)量變化，提出了保證蒙特卡羅輸運計算精度的密度修正。

除UMV修正和密度修正外，選擇合適的截面庫也很重要，包括輸運計算的連續(xù)截面庫和燃耗計算中采用的燃耗截面庫。在MOI研制過程中，對截面庫也進(jìn)行過一些改進(jìn)，如通過修正燃耗截面庫中Nd的裂變產(chǎn)額，改善143Nd的核密度。

圖8示出IFBA例題誤差修正后的計算結(jié)果。

圖8 誤差修正后kinf隨燃耗的變化

2 程序驗證

驗證基準(zhǔn)題[11]采用的柵元模型來自西屋壓水堆17×17組件，如圖9所示。采用ThO2-UO2芯塊來替代UO2燃料芯塊，實際密度是理論密度的94%；重核素中Th的質(zhì)量占75%，U的質(zhì)量占25%(235U的富集度為19.5%)。表1、2列出西屋壓水堆組件單柵元模型的詳細(xì)參數(shù)，本工作的基準(zhǔn)計算采用熱態(tài)滿功率參數(shù)。

圖9 單柵元模型

表1 單柵元模型參數(shù)

表2 熱態(tài)滿功率下柵元的初始成分

采用基準(zhǔn)題計算的棒柵元kinf隨燃耗的變化示于圖10。作為比較，圖10中同時給出CASMO-4、MIT和INEEL采用MOCUP的計算結(jié)果。由圖10可見，MOI計算的kinf與MIT-MOCUP的結(jié)果幾乎完全一致。不同程序計算的棒柵元各核素含量與CASMO-4的相對誤差示于圖11。由圖11可見，MOI計算的核密度與CASMO-4的相對誤差均在7%以內(nèi)，而MOCUP計算的個別核密度與CASMO-4的相對誤差超過10%。

由圖10、11可知，誤差修正后，MOI的計算精度得到了改進(jìn)，單柵元模型計算的有效增殖因數(shù)達(dá)到了國際同類程序的水平，而核密度則更接近成熟的商業(yè)程序CASMO-4，為MOI的實際應(yīng)用奠定了基礎(chǔ)。

圖10 基準(zhǔn)題棒柵元kinf隨燃耗的變化

圖11 不同程序計算的棒柵元各核素含量與CASMO-4計算結(jié)果的相對誤差

3 結(jié)論

本工作對影響蒙卡燃耗計算誤差的原因進(jìn)行了分類，對蒙特卡羅方法統(tǒng)計誤差的傳遞效應(yīng)進(jìn)行了研究。以MOI為基礎(chǔ)，采用UMV修正、密度修正等措施，在一定程度上降低了蒙卡燃耗計算中的誤差，有利于蒙卡燃耗計算程序在反應(yīng)堆工程設(shè)計驗證中的應(yīng)用。

參考文獻(xiàn)：

[1] BAKKARI B E, ELBARDOUNI T, NACIR B, et al. Accuracy assessment of a new Monte Carlo based burnup computer code[J]. Annals of Nuclear Energy, 2012, 45: 29-36.

[2] TAKEDA T, HIROKAWA N, NODA T. Estimation of error propagation in Monte-Carlo burnup calculations[J]. J Nucl Sci Tech, 1999, 36: 738-745.

[3] TOHJOH M, ENDO T, WATANABE M, et al. Effect of error propagation of nuclide number densities on Monte Carlo burn-up calculations[J]. Annals of Nuclear Energy, 2006, 33: 1 424-1 436.

[4] DUMONTEIL E, DIOP C M. Biases and statistical errors in Monte Carlo burnup calculations: An unbiased stochastic scheme to solve Boltzmann/Bateman coupled equations[J]. Nuclear Science and Engineering, 2011, 167: 165-170.

[5] YAMAMOTO A, TATSUMI M, SUGIMURA N. Numerical solution of stiff burnup equation with short half lived nuclides by the Krylov subspace method[J]. Journal of Nuclear Science and Technology, 2007, 44(2): 147-154.

[6] ISOTALO A E, AARNIO P A. Substep methods for burnup calculations with Bateman solutions[J]. Annals of Nuclear Energy, 2011, 38: 2 509-2 514.

[7] ISOTALO A E, AARNIO P A. Higher order methods for burnup calculations with Bateman solutions[J]. Annals of Nuclear Energy, 2011, 38: 1 987-1 995.

[8] XU Zhiwen, HEJZLAR P, DRISCOLL M J, et al. An improved MCNP-Origen depletion program (MCODE) and its verification for high-burnup applications[C]∥PHYSOR 2002. Seoul: [s. n.], 2002.

[9] 余綱林，王侃. MCNP和ORIGEN2 耦合系統(tǒng)(MCBurn)的研究[D]. 北京：清華大學(xué)，2002.

[10] 佘頂，王侃，余綱林. 堆用蒙卡程序燃耗計算功能開發(fā)[J]. 核動力工程，2012，33(3)：1-5.

SHE Ding, WANG Kan, YU Ganglin. Development of burnup calculation function in reactor Monte Carlo code RMC[J]. Nuclear Power Engineering, 2012, 33(3): 1-5(in Chinese).

[11] WEAVER K D, ZHAO X, PILAT E E, et al. A PWR thorium pin cell burnup benchmark[C]∥PHYSOR 2002. Seoul: [s. n.], 2002.