基于RMC程序的keff對(duì)核數(shù)據(jù)的敏感性分析

丘意書(shū)，余健開(kāi)，梁金剛，王 侃

(清華大學(xué) 工程物理系，北京 100084)

基于RMC程序的keff對(duì)核數(shù)據(jù)的敏感性分析

丘意書(shū)，余健開(kāi)，梁金剛，王 侃

(清華大學(xué) 工程物理系，北京 100084)

適用于連續(xù)能量蒙特卡羅程序的敏感性分析方法是當(dāng)前的研究熱點(diǎn)。本文建立了5種不同反應(yīng)類型的敏感性系數(shù)的計(jì)算公式，對(duì)當(dāng)前應(yīng)用廣泛的反復(fù)裂變幾率法的理論基礎(chǔ)及算法進(jìn)行了分析。分別使用RMC程序和MCNP6程序計(jì)算了keff對(duì)核數(shù)據(jù)的敏感性系數(shù)，計(jì)算結(jié)果吻合良好。本文結(jié)果表明RMC程序初步具備了敏感性分析的功能。

敏感性分析；反復(fù)裂變幾率法；連續(xù)能量；RMC程序

當(dāng)使用測(cè)量?jī)x器或數(shù)值方法得到某個(gè)物理量的估計(jì)值時(shí)，需要知道它離真值多近，估計(jì)值和真值的差異用誤差來(lái)表示。然而，真值往往是未知的或不可知的，一般只能估計(jì)誤差的大小，對(duì)誤差的估計(jì)稱之為不確定度。不確定性分析的任務(wù)就是分析由輸入(例如核數(shù)據(jù)因測(cè)量引起)的不確定度導(dǎo)致計(jì)算結(jié)果(keff等)的不確定度。

不確定性分析的方法一般包括兩類。第1類是隨機(jī)抽樣法，根據(jù)方差及協(xié)方差的大小抽樣輸入?yún)?shù)，產(chǎn)生無(wú)窮多個(gè)輸入系列，進(jìn)行輸運(yùn)計(jì)算，然后對(duì)計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)分析，得到keff的不確定度，代表程序?yàn)榈聡?guó)核設(shè)備與反應(yīng)堆安全研究協(xié)會(huì)(GRS)研制的XSUSA程序[1]。該方法簡(jiǎn)單，容易實(shí)現(xiàn)，但需非常大的計(jì)算量。第2類方法是基于敏感性分析的方法。該方法需先計(jì)算出keff對(duì)核截面的敏感性系數(shù)，再結(jié)合截面的協(xié)方差數(shù)據(jù)庫(kù)計(jì)算出keff的不確定度。這類方法計(jì)算效率高，是當(dāng)前用于keff對(duì)核數(shù)據(jù)不確定性分析的主要方法，代表程序?yàn)槊绹?guó)橡樹(shù)嶺國(guó)家實(shí)驗(yàn)室(ORNL)開(kāi)發(fā)的SCALE6.1程序包[2]中一維敏感性與不確定性分析模塊TSUNAMI-1D和三維敏感性與不確定性分析模塊TSUNAMI-3D。TSUNAMI-1D和TSUNAMI-3D的輸運(yùn)計(jì)算分別使用的是SN程序XSDRNPM和蒙特卡羅程序KENO，它們的共同點(diǎn)是基于多群截面數(shù)據(jù)庫(kù)，需要執(zhí)行并群和共振自屏計(jì)算，因而必須計(jì)算隱式敏感性系數(shù)，以考慮多群近似對(duì)敏感性系數(shù)的影響。另外，需要執(zhí)行一次前向計(jì)算和伴隨計(jì)算分別獲得通量和伴隨通量的信息。

為使計(jì)算流程更加簡(jiǎn)單及使計(jì)算結(jié)果更加可靠，使用連續(xù)能量蒙特卡羅程序計(jì)算keff對(duì)核數(shù)據(jù)的敏感性系數(shù)，成為近幾年蒙特卡羅方法研究的熱點(diǎn)。然而，現(xiàn)有基于多群蒙特卡羅程序的敏感性分析方法無(wú)法適用于連續(xù)能量蒙特卡羅程序中，因?yàn)樵趫?zhí)行伴隨蒙特卡羅計(jì)算時(shí)對(duì)連續(xù)能量的中子散射算符進(jìn)行轉(zhuǎn)置非常困難。因此，需要提出新的適用于連續(xù)能量蒙特卡羅框架的計(jì)算敏感性系數(shù)的方法。這些方法包括微分算符抽樣法(DOS)[3]、反復(fù)裂變幾率法(IFP)[4]、次級(jí)粒子貢獻(xiàn)法(Contributon)[5]、CLUTCH法[6]及Contributon-IFP混合法[7]。其中，反復(fù)裂變幾率法最先在計(jì)算動(dòng)態(tài)參數(shù)上獲得應(yīng)用[8]，被用于敏感性分析研究[6，9-10]。由于該方法物理概念清晰，計(jì)算精度高，幾乎被現(xiàn)在所有具備敏感性分析功能的連續(xù)能量蒙特卡羅程序(MCNP6.1[11]、SCALE[7]及McCARD[12]程序)所采用。因而，本文選用該方法進(jìn)行keff敏感性分析。

1 理論基礎(chǔ)

基于一階線性微擾理論[13]，從中子輸運(yùn)方程和伴隨方程可推導(dǎo)得到由核數(shù)據(jù)的擾動(dòng)引起的keff的擾動(dòng)為：

(1)

其中：Ψ為通量；Ψ*為伴隨通量；k表示keff；Σt為宏觀總截面；S為散射算符；F為裂變算符；尖括號(hào)表示對(duì)相空間進(jìn)行積分。

敏感性系數(shù)定義為某個(gè)核數(shù)據(jù)的相對(duì)變化量所引起的keff的相對(duì)變化量，如式(2)所示。

(2)

根據(jù)敏感性系數(shù)的定義及式(1)可得：

(3)

(4)

下面分別給出上述5種不同反應(yīng)類型的敏感性系數(shù)的計(jì)算公式。

keff對(duì)吸收反應(yīng)的敏感性系數(shù)為：

(5)

(6)

其中：g為能群號(hào)；z為所關(guān)心的區(qū)域；V為體積；Σa為吸收截面；Σf為裂變截面；E′為引起裂變的中子的能量；Ω′為引起裂變的中子的方向角；E為裂變產(chǎn)生的中子的能量；Ω為裂變產(chǎn)生的中子的方向角；D為伴隨通量加權(quán)的裂變?cè)错?xiàng)。

keff對(duì)散射反應(yīng)的敏感性系數(shù)為：

(7)

其中，Σs，g為能群g的散射截面。式(7)中第1項(xiàng)代表了玻爾茲曼方程中的散射項(xiàng)對(duì)所求的敏感性系數(shù)的貢獻(xiàn)，第2項(xiàng)代表了玻爾茲曼方程中的碰撞項(xiàng)對(duì)該敏感性系數(shù)的貢獻(xiàn)。

keff對(duì)裂變反應(yīng)的敏感性系數(shù)為：

(8)

dΩ′dΩdE′dVdE

(9)

keff對(duì)裂變中子能譜χ的敏感性系數(shù)為：

(10)

與前4類敏感性系數(shù)不同，計(jì)算keff對(duì)χ的敏感性系數(shù)需對(duì)出射能量劃分網(wǎng)格。

2 反復(fù)裂變幾率法

2.1 理論基礎(chǔ)

早在1964年，Hurwitz[14]就已指出，反復(fù)裂變幾率與伴隨通量呈正比，其物理含義為引入反應(yīng)堆的某個(gè)中子在無(wú)窮遠(yuǎn)代產(chǎn)生的裂變次數(shù)。

伴隨方程可寫為：

(11)

其中，L*為伴隨算符。

使用源迭代法求解伴隨通量，第n代的伴隨通量的解為：

(12)

(13)

(14)

其中，Ln(P→P′)表示相空間P的父代中子在P′處產(chǎn)生的第n代裂變中子。

因此，伴隨通量可用反復(fù)裂變幾率來(lái)表示。經(jīng)數(shù)值驗(yàn)證，在連續(xù)能量蒙特卡羅程序中，一般取n=10可使大部分系統(tǒng)的反復(fù)裂變幾率收斂。

2.2 算法實(shí)現(xiàn)

反復(fù)裂變幾率法是使用反復(fù)裂變幾率作為伴隨通量的估計(jì)值的方法，該方法直接在前向輸運(yùn)計(jì)算中統(tǒng)計(jì)得到反復(fù)裂變幾率的大小，不需求解共軛方程而額外執(zhí)行1次伴隨計(jì)算。根據(jù)反復(fù)裂變幾率的物理含義，可人為地將蒙特卡羅模擬的活躍代劃分為連續(xù)的塊，每個(gè)塊的大小為n，應(yīng)使得伴隨通量收斂。每個(gè)塊的第1代稱為初始代，初始代的裂變中子為父代中子，這些父代中子會(huì)產(chǎn)生子代中子，足夠多的過(guò)渡代后(共n-2代)，父代中子產(chǎn)生的子代中子的數(shù)量會(huì)趨于穩(wěn)定，此時(shí)，對(duì)應(yīng)的活躍代稱為漸進(jìn)代，也就是每塊的最后1代。根據(jù)敏感性系數(shù)的計(jì)算公式，在初始代統(tǒng)計(jì)特定核素、特定反應(yīng)類型的反應(yīng)率，在漸進(jìn)代統(tǒng)計(jì)對(duì)應(yīng)的反復(fù)裂變幾率。具體的實(shí)現(xiàn)可使用伴隨權(quán)重計(jì)數(shù)法，不需直接計(jì)算出伴隨通量，因而無(wú)須對(duì)空間、能群、角度進(jìn)行離散，計(jì)算精確高。

伴隨權(quán)重計(jì)數(shù)原理如圖1所示。在初始代，需定義父代中子的軌跡，它定義為父代中子出生地點(diǎn)和源中子之間的軌跡，如圖1所示，源中子發(fā)生了兩次裂變，因而可定義兩條父代中子軌跡，分別是T1及T1+T2。根據(jù)所統(tǒng)計(jì)敏感性系數(shù)的類型，可分別獲得δT1和δ(T1+T2)的貢獻(xiàn)，再結(jié)合徑跡長(zhǎng)度估計(jì)法可獲得初始代中反應(yīng)率的估計(jì)值。其中，如果發(fā)生了所統(tǒng)計(jì)的敏感性系數(shù)的反應(yīng)，δ等于1，否則為零。同時(shí)，為了跟蹤后代中子，對(duì)每個(gè)父代中子進(jìn)行編號(hào)。圖1中兩個(gè)祖先中子的編號(hào)分別為1和2。在過(guò)渡代，通過(guò)讓所有后代中子繼續(xù)父代中子編號(hào)的方法，跟蹤父代中子產(chǎn)生的子代中子。在漸進(jìn)代，統(tǒng)計(jì)這些父代中子所對(duì)應(yīng)的子代中子(具有相同編號(hào))引起的裂變以作為伴隨權(quán)重的大小，圖1中分別對(duì)應(yīng)Q1及Q2+Q3。則兩個(gè)父代中子的計(jì)數(shù)分別為T1Q1和(T1+T2)、(Q2+Q3)。

圖1 伴隨權(quán)重計(jì)數(shù)原理Fig.1 Scheme of adjoint weighted tally

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

根據(jù)上述基本原理，基于清華大學(xué)工程物理系研制的自主堆用連續(xù)能量蒙特卡羅程序RMC[15]，使用反復(fù)裂變幾率法開(kāi)發(fā)了keff對(duì)核數(shù)據(jù)的敏感性分析功能，并將RMC的計(jì)算結(jié)果與MCNP6的進(jìn)行比較。RMC的版本為RMC2.0，MCNP6為2013年發(fā)布的MCNP6.1。選取一個(gè)3×3燃料柵格作為驗(yàn)證基準(zhǔn)題[16]，它的幾何構(gòu)造如圖2所示。使用了27個(gè)不同的核素對(duì)該基準(zhǔn)題進(jìn)行建模。表1列出RMC和MCNP6的keff的計(jì)算結(jié)果，計(jì)算條件為每代200 000個(gè)粒子，非活躍代275代，一共20 275代，RMC和MCNP6的keff小于1 pcm。

為了統(tǒng)計(jì)敏感性系數(shù)，將活躍代劃分成塊，塊大小取為10。表2列出RMC和MCNP6計(jì)算得到的重要核數(shù)據(jù)的敏感性系數(shù)。表2中，C/E定義為RMC與MCNP6的敏感性系數(shù)的比值。從表2可見(jiàn)，C/E的比值整體在1附近，兩個(gè)程序計(jì)算結(jié)果總體吻合良好。keff對(duì)235U的平均裂變中子數(shù)最為敏感，約0.93。另外，RMC和MCNP6均滿足keff對(duì)所有裂變核素的平均裂變中子數(shù)的敏感性系數(shù)之和等于1的檢驗(yàn)標(biāo)準(zhǔn)。

圖2 基準(zhǔn)題幾何Fig.2 Geometry of benchmark

表1 keff的計(jì)算結(jié)果Table 1 Calculation results of keff

表2中，兩個(gè)程序所得的敏感性系數(shù)的相對(duì)偏差相當(dāng)，表明它們的計(jì)數(shù)效率相當(dāng)。與一般的蒙特卡羅計(jì)數(shù)不同，使用蒙特卡羅統(tǒng)計(jì)得到的敏感性系數(shù)并非其值越小，其相對(duì)偏差就一定越大。相對(duì)偏差的大小既與敏感性系數(shù)本身的大小有關(guān)，同時(shí)也與計(jì)數(shù)的方法有關(guān)。對(duì)于吸收反應(yīng)類型，由于其敏感性系數(shù)只由碰撞項(xiàng)構(gòu)成，兩個(gè)程序均采用徑跡長(zhǎng)度法，抽樣效率較高。對(duì)于散射及裂變反應(yīng)類型，其敏感性系數(shù)不僅包括碰撞項(xiàng)，還包括散射或裂變項(xiàng)，RMC使用碰撞估計(jì)法對(duì)它們進(jìn)行統(tǒng)計(jì)，而MCNP6采用平均估計(jì)法統(tǒng)計(jì)。由于采用不同的統(tǒng)計(jì)方法，敏感性系數(shù)的相對(duì)偏差也就不同。例如，盡管表2中keff對(duì)1H彈性散射截面的敏感性系數(shù)比(n,γ)的敏感性系數(shù)要大，但是其相對(duì)偏差也較大。

表2 不同核數(shù)據(jù)的敏感性系數(shù)Table 2 Sensitivity coefficients for different types of nuclear data

表3列出RMC和MCNP6計(jì)算得到的重要核素的總截面的敏感性系數(shù)。由表3可見(jiàn)，keff對(duì)235U、238U、1H 和10B 4種核素最為敏感。其中，keff對(duì)235U和1H的總截面的敏感性系數(shù)是正值，而對(duì)238U和10B的總截面的敏感性系數(shù)是負(fù)值。因?yàn)?35U所有截面中對(duì)keff影響最大的是裂變截面，1H所有截面中對(duì)keff影響最大的是散射截面，它們均與keff正相關(guān)，而238U和10B所有截面中對(duì)keff影響最大的是(n,γ)截面，它與keff負(fù)相關(guān)。

圖3示出235U和238U的重要敏感性系數(shù)曲線，能量網(wǎng)格的劃分與SCALE6程序中的AMPX多群數(shù)據(jù)庫(kù)中的中子238群能群結(jié)構(gòu)一致。從圖3a可看出，keff對(duì)235U裂變截面的敏感性系數(shù)與它對(duì)235U平均裂變中子數(shù)的趨勢(shì)基本一致，在低能區(qū)出現(xiàn)峰值，這與235U容易發(fā)生裂變的能量范圍一致，而(n,γ)截面的敏感性系數(shù)在熱能區(qū)出現(xiàn)谷值，表明在低能區(qū)235U的吸收截面較大。從圖3b可看出，keff對(duì)238U裂變截面的敏感性系數(shù)與它對(duì)238U平均裂變中子數(shù)的趨勢(shì)基本一致，在快中子區(qū)出現(xiàn)峰值，這與238U容易發(fā)生裂變的能量范圍一致，(n,γ)截面的敏感性系數(shù)在共振區(qū)有強(qiáng)烈的波動(dòng)，反映了238U的共振峰對(duì)keff的影響。RMC的計(jì)算結(jié)果和MCNP6的吻合良好。

表3 總截面的敏感性系數(shù)Table 3 Sensitivity coefficients for total cross sections

圖3 235U(a)和238U(b)的重要敏感性系數(shù)曲線Fig.3 Curves of some important sensitivity coefficients of 235U (a) and 238U (b)

4 結(jié)論與展望

本文討論了敏感性系數(shù)的理論基礎(chǔ)及反復(fù)裂變幾率法的基本原理，基于自主堆用蒙特卡羅程序RMC，開(kāi)發(fā)了連續(xù)能量蒙特卡羅計(jì)算中keff對(duì)核數(shù)據(jù)的敏感性分析功能。選用了3×3燃料柵格臨界基準(zhǔn)題，對(duì)RMC及MCNP6的計(jì)算結(jié)果進(jìn)行了比較。RMC和MCNP6的keff偏差小于1 pcm，RMC計(jì)算得到的敏感性系數(shù)與MCNP6計(jì)算的比值基本在1左右，吻合良好。敏感性分析是不確定性分析的基礎(chǔ)，本文工作為后續(xù)的不確定性分析打下了基礎(chǔ)。下一步工作將結(jié)合協(xié)方差數(shù)據(jù)庫(kù)，使用RMC進(jìn)行keff對(duì)核數(shù)據(jù)的不確定性分析。

[1] ZWERMANN W, GALLNER L, KLEIN M, et al. XSUSA solution for the PWR Pincell burnup benchmark[C]∥Proceedings of the 6th Workshop on OECD Benchmark for Uncertainty Analysis in Best-Estimate Modeling for Design, Operation and Safety Analysis of LWRs (UAM-6). Karlsruhe, Germany: [s. n.], 2012.

[2] SCALE: A comprehensive modeling and simulation suite for nuclear safety analysis and design, Version 6.1, ORNL/TM-2005/39[R]. USA: Oak Ridge National Laboratory, 2011.

[3] KIEDROWSKI B C, BROWN F B. Comparison of the Monte Carlo adjoint-weighted and differential operator perturbation methods[R]. USA: Los Alamos National Laboratory, 2010.

[4] NAUCHI Y, KAMEYAMA T. Development of calculation technique for iterated fission probability and reactor kinetic parameters using continuous-energy Monte Carlo method[J]. J Nucl Sci Technol, 2010, 47(11): 977-990.

[5] REARDEN B T, PERFETTI C M, WILLIAMS M L, et al. SCALE sensitivity calculations using contributon theory[C]∥ Proceedings of Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo 2010 (SNA+MC2010). Tokyo, Japan: [s. n.], 2010.

[6] PERFETTI C M, REARDEN B T. Development of a SCALE tool for continuous-energy eigenvalue sensitivity coefficient calculations[C]∥Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo 2013 (SNA+MC 2013). Paris, France: [s. n.], 2013.

[7] PERFETTI C M. Advanced Monte Carlo methods for continuous-energy eigenvalue sensitivity coefficient calculation[D]. USA: University of Michigan, 2012.

[8] KIEDROWSKI B C, BROWN F B, WILSON P P H. Adjoint-weighted tallies fork-eigenvalue calculations with continuous-energy Monte Carlo[J]. Nuclear Science and Engineering, 2011, 168(3): 226-241.

[9] KIEDROWSKI B C, BROWN F B. Adjoint-basedk-eigenvalue sensitivity coefficients to nuclear data using continuous-energy Monte Carlo[J]. Nuclear Science and Engineering, 2012, 174: 227-244.

[10]SHIM H J, KIM C H. Adjoint sensitivity and uncertainty analyses in Monte Carlo forward calculations[J]. Journal of Nuclear Science and Technology, 2011, 5: 1 453-1 461.

[11]KIEDROWSKI B C. MCNP6.1k-eigenvalue sensitivity capability: A users guide, LA-UR-13-22251[R]. USA: LANL, 2013.

[12]SHIM H J, HAN B S, JUNG J S, et al. MCCARD: Monte Carlo code for advanced reactor design and analysis[J]. Nuclear Engineering and Technology, 2012, 44(2): 161-176.

[13]REARDEN B T. Perturbation theory eigenvalue sensitivity analysis with Monte Carlo techniques[J]. Nucl Sci Eng, 2004, 146: 367-382.

[14]HURWITZ H. Physical interpretation of the adjoint flux: iterated fission probability[M]∥Naval reactor physics handbook, Vol. Ⅰ. RADKOWSKY A. USA: Atomic Energy Commission, 1964: 864-869.

[15]WANG Kan, LI Zeguang, SHE Ding, et al. RMC: A Monte Carlo code for reactor core analysis[C]∥Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo 2013 (SNA+MC 2013). Paris, France: [s. n.], 2013.

[16]DALLE H M. Monte Carlo burnup simulation of the Takahama-3 benchmark experiment[C]∥2009 International Nuclear Atlantic Conference: INAC 2009. Brazil: [s. n.], 2009.

keffSensitivity Analysis to Nuclear Data with RMC Code

QIU Yi-shu, YU Jian-kai, LIANG Jin-gang, WANG Kan

(DepartmentofEngineeringPhysics,TsinghuaUniversity,Beijing100084,China)

Methods suitable for sensitivity analysis in continuous-energy Monte Carlo codes become a research hotspot in the field of reactor physics. In this work, the formulas of sensitivity coefficients of five different reaction types were established. Then, the theoretical basis and the algorithm of the iterated fission probability method which was used widely currently were discussed. Furthermore, two Monte Carlo codes, RMC and MCNP6, were used to compute eigenvalue sensitivity coefficients to nuclear data. The agreement between RMC and MCNP6 is well. The results indicate that RMC is capable to perform sensitivity analysis preliminarily.

sensitivity analysis; iterated fission probability method; continuous-energy; RMC code

2014-06-18;

2014-12-02

國(guó)家自然科學(xué)基金資助項(xiàng)目(11475098)

丘意書(shū)(1990—)，男，廣東河源人，博士研究生，反應(yīng)堆物理專業(yè)

TL32

A

1000-6931(2015)10-1821-07

10.7538/yzk.2015.49.10.1821