黃金龍，曹良志,*，賀清明，李 捷，沈 煒,2

(1.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049；2.加拿大CANDU Owners Group)

核數(shù)據(jù)的敏感性分析是核數(shù)據(jù)不確定度量化的基礎(chǔ)。相比于基于確定論程序的多群核數(shù)據(jù)敏感性分析[1]，基于蒙特卡羅程序的連續(xù)能量核數(shù)據(jù)敏感性分析具有幾何適應(yīng)性強(qiáng)、計(jì)算精度高等優(yōu)勢(shì)。近10年中，連續(xù)能量核數(shù)據(jù)敏感性分析在國(guó)際上取得了巨大進(jìn)展。國(guó)際主流的蒙特卡羅程序，如MCNP[2]、Serpent[3]、SCALE[4]、McCARD[5]、OpenMC[6]等已具備連續(xù)能量核數(shù)據(jù)敏感性分析的功能，國(guó)內(nèi)的蒙特卡羅程序RMC[7]也已具備此功能。從方法上看，MCNP首先采用反復(fù)裂變幾率(IFP)方法進(jìn)行連續(xù)能量核數(shù)據(jù)敏感性分析，此方法的優(yōu)點(diǎn)是計(jì)算精度高，缺點(diǎn)是占用內(nèi)存大、統(tǒng)計(jì)漲落較大。SCALE在此基礎(chǔ)上提出了CLUTCH方法，有效降低了存儲(chǔ)內(nèi)存，但計(jì)算過(guò)程需要?jiǎng)澐志W(wǎng)格，計(jì)算存在近似。RMC提出了超歷史法[8]以及基于裂變的CLUTCH等方法[9]提高計(jì)數(shù)效率。Serpent提出了碰撞歷史法[3]，此方法是一種極具創(chuàng)新性的方法，可用于求解特征值靈敏度系數(shù)、廣義靈敏度系數(shù)、共振參數(shù)的靈敏度系數(shù)和中子代時(shí)間的靈敏度系數(shù)等，同時(shí)Serpent還采用了擴(kuò)展廣義微擾理論[10]進(jìn)行靈敏度系數(shù)的計(jì)算，該方法基于基函數(shù)展開，得到連續(xù)的函數(shù)形式的靈敏度系數(shù)公式，并基于連續(xù)形式的協(xié)方差信息，得到連續(xù)形式的響應(yīng)參數(shù)的不確定度，使得靈敏度系數(shù)和不確定度結(jié)果實(shí)現(xiàn)“連續(xù)”。

本文在西安交通大學(xué)核工程計(jì)算物理實(shí)驗(yàn)室自主開發(fā)的蒙特卡羅程序NECP-MCX[11]中開發(fā)連續(xù)能量核數(shù)據(jù)敏感性分析的功能模塊。采用反復(fù)裂變幾率方法求解靈敏度系數(shù)，并使用稀疏矩陣的存儲(chǔ)方式以及重疊塊法分別降低存儲(chǔ)內(nèi)存和提高計(jì)數(shù)效率，最后開展AP1000堆芯keff對(duì)連續(xù)能量核數(shù)據(jù)的敏感性分析研究。

1 理論方法

本文從中子輸運(yùn)方程出發(fā)，推導(dǎo)keff對(duì)連續(xù)能量核數(shù)據(jù)的靈敏度系數(shù)計(jì)算公式，并結(jié)合蒙特卡羅程序模擬的特點(diǎn)，利用反復(fù)裂變幾率方法得到共軛通量，從而在前向蒙特卡羅模擬中得到keff對(duì)連續(xù)能量核數(shù)據(jù)的靈敏度系數(shù)。

1.1 一階微擾理論

穩(wěn)態(tài)中子輸運(yùn)方程可表示為：

(1)

其中：B為輸運(yùn)-碰撞項(xiàng)算子；ψ為中子通量密度；k為有效增殖因子；F為產(chǎn)生項(xiàng)算子。

輸運(yùn)-碰撞項(xiàng)算子B表示為：

Ω′→Ω)ψ(r,Ω′,E′)dE′dΩ′

(2)

產(chǎn)生項(xiàng)算子F表示為：

ψ(r,Ω′,E′)dE′dΩ′

(3)

同樣，穩(wěn)態(tài)共軛中子輸運(yùn)方程可表示為：

(4)

共軛輸運(yùn)-碰撞項(xiàng)算子B*表示為：

Ω→Ω′)ψ*(r,Ω′,E′)dE′dΩ′

(5)

共軛產(chǎn)生項(xiàng)算子F*表示為：

ψ*(r,Ω′,E′)dE′dΩ′

(6)

引入擾動(dòng)：

F′=F+δF

B′=B+δB

k′=k+δk

ψ′=ψ+δψ

(7)

則擾動(dòng)后的中子輸運(yùn)方程表示為：

(8)

將式(8)兩邊與ψ*在相空間中內(nèi)積，得到：

〈ψ*,k′B′ψ′〉=〈ψ*,F′ψ′〉

(9)

將式(1)兩邊與ψ′在相空間中內(nèi)積，得到：

〈ψ*,(kB-F)ψ′〉=〈ψ′,(kB*-F*)ψ*〉

(10)

根據(jù)式(9)、(10)以及共軛算子的性質(zhì)，采用一階近似可得到：

(11)

靈敏度系數(shù)的定義為：

(12)

將式(11)代入式(12)，可得到：

(13)

對(duì)式(13)分子中的3項(xiàng)以及分母分別展開，可得到裂變項(xiàng)：

(14)

散射項(xiàng)：

Ω′→Ω)ψ(r,Ω′,E′)dΩ′dΩdr

(15)

碰撞項(xiàng)：

Σx(r,E)ψ(r,Ω,E)dΩdr

(16)

分母表示為：

Σf(r,Ω′,E′)ψ(r,Ω′,E′)dE′dΩ′·

(17)

針對(duì)具體的核數(shù)據(jù)，根據(jù)式(14)～(17)的組合可得到keff對(duì)該核數(shù)據(jù)的靈敏度系數(shù)。

1.2 反復(fù)裂變幾率

從式(14)～(17)可看出，求解連續(xù)能量核數(shù)據(jù)的靈敏度系數(shù)需得到共軛通量以及核反應(yīng)率，統(tǒng)計(jì)核反應(yīng)率在蒙特卡羅程序中易實(shí)現(xiàn)，因此求解靈敏度系數(shù)的關(guān)鍵在于求解共軛通量。

Hurwitz[12]首先解釋了反復(fù)裂變幾率的物理含義：在相空間(r，Ω，E)處引入1個(gè)中子，稱為祖先中子，經(jīng)過(guò)數(shù)代蒙特卡羅模擬后，由該祖先中子產(chǎn)生的后代中子數(shù)目會(huì)達(dá)到穩(wěn)定，后代中子的數(shù)目即為(r，Ω，E)處的反復(fù)裂變幾率，具體過(guò)程如圖1所示?？勺C明，相空間中某點(diǎn)(r，Ω，E)的反復(fù)裂變幾率正比于該點(diǎn)的共軛通量，由于式(13)中分子和分母同時(shí)包含共軛通量，因此可直接用反復(fù)裂變幾率代替式(13)中的共軛通量。基于反復(fù)裂變幾率方法可直接在前向蒙特卡羅模擬中得到共軛通量。蒙特卡羅模擬中使祖先中子產(chǎn)生的后代中子數(shù)達(dá)到平衡所需要的蒙特卡羅代數(shù)稱為塊，塊的大小由問(wèn)題決定，一般塊的大小取5～20。在1個(gè)塊中，將蒙特卡羅代分為初始代、過(guò)渡代和漸近代，在初始代中統(tǒng)計(jì)核反應(yīng)率，并對(duì)祖先中子進(jìn)行編號(hào)，過(guò)渡代中傳遞祖先中子的編號(hào)信息，漸近代統(tǒng)計(jì)后代中子的數(shù)目。因此，1個(gè)塊由數(shù)個(gè)蒙特卡羅代組成，經(jīng)過(guò)1個(gè)塊的模擬，可得到1次靈敏度系數(shù)的計(jì)數(shù)。

圖1 反復(fù)裂變幾率原理示意圖Fig.1 Schematic diagram of iterated fission probability

由于需要統(tǒng)計(jì)數(shù)代后的中子數(shù)目，需在數(shù)代之后才能得到共軛通量的計(jì)數(shù)，無(wú)法在蒙特卡羅模擬的當(dāng)代得到共軛通量的計(jì)數(shù)，且核反應(yīng)率與共軛通量的乘積與粒子相關(guān)，因此需存儲(chǔ)與模擬粒子數(shù)相關(guān)的核反應(yīng)率計(jì)數(shù)數(shù)組。當(dāng)粒子數(shù)較大、計(jì)算的靈敏度系數(shù)較多時(shí)，內(nèi)存占用非常大。對(duì)于VERA 1B算例，總代數(shù)為500代，其中非活躍代設(shè)置為200代。數(shù)值結(jié)果表明，當(dāng)粒子數(shù)大于100時(shí)，核反應(yīng)率數(shù)組非零值比例達(dá)到穩(wěn)定，即不再變化，因此本測(cè)試粒子數(shù)取為100，在2群和44群能群劃分下統(tǒng)計(jì)不同核素不同反應(yīng)道的核反應(yīng)率，其非零值比例列于表1。

表1 核反應(yīng)率計(jì)數(shù)數(shù)組非零值比例Table 1 None-zero value rate in reaction rate array

從表1可看出，核反應(yīng)率計(jì)數(shù)數(shù)組中非零值的比例與能群結(jié)構(gòu)有關(guān)。在2群能群結(jié)構(gòu)下，由于能群劃分較粗，非零值比例較大，達(dá)到82%。而在44群能群結(jié)構(gòu)下，核反應(yīng)率計(jì)數(shù)數(shù)組中非零值的比例最大僅為15%，因此核反應(yīng)率數(shù)組中存儲(chǔ)的大部分?jǐn)?shù)據(jù)是對(duì)結(jié)果沒(méi)有明顯意義的零值，可采用稀疏矩陣的方式存儲(chǔ)核反應(yīng)率計(jì)數(shù)數(shù)組，以達(dá)到降低內(nèi)存的效果。

在反復(fù)裂變幾率方法中，不同的塊依次相連稱為等待法，由于需用漸近代統(tǒng)計(jì)得到的共軛通量對(duì)初始代統(tǒng)計(jì)得到的核反應(yīng)率進(jìn)行加權(quán)，才能得到1次靈敏度系數(shù)的計(jì)數(shù)，因此等待法計(jì)數(shù)效率較低。而重疊塊法[13]中塊與塊之間有重疊，第i代的中子，既處于第i塊的初始代，也是第i-B+2塊至第i-1塊的過(guò)渡代(B表示塊的大小，i>B)，同時(shí)也是第i-B+1塊的漸近代，塊的數(shù)量明顯增多，因此能提高計(jì)數(shù)效率。等待法和重疊塊法示意圖如圖2所示。如果將蒙特卡羅模擬的條件設(shè)置為300代活躍代、塊的大小取10，采用等待法進(jìn)行計(jì)算，靈敏度系數(shù)計(jì)數(shù)次數(shù)為30，而采用重疊塊法進(jìn)行計(jì)算，靈敏度系數(shù)的計(jì)數(shù)次數(shù)為291。

圖2 等待法和重疊塊法區(qū)別Fig.2 Difference of waiting method and overlapping block method

由于蒙特卡羅模擬中隨機(jī)變量的統(tǒng)計(jì)標(biāo)準(zhǔn)差σ滿足：

(18)

其中，N為試驗(yàn)次數(shù)。因此，采用重疊塊法計(jì)數(shù)的試驗(yàn)次數(shù)N大于等待法，重疊塊法計(jì)數(shù)的統(tǒng)計(jì)漲落小于等待法。

2 程序驗(yàn)證

根據(jù)前文介紹的基于重疊塊法的反復(fù)裂變幾率方法、稀疏矩陣存儲(chǔ)方式，在蒙特卡羅程序NECP-MCX上完成keff對(duì)連續(xù)能量核數(shù)據(jù)敏感性分析功能模塊的開發(fā)，對(duì)VERA 2B算例進(jìn)行驗(yàn)證，并將NECP-MCX的靈敏度系數(shù)計(jì)算結(jié)果與MCNP6進(jìn)行對(duì)比。

2.1 算例設(shè)置

本文選取的驗(yàn)證算例取自VERA系列基準(zhǔn)題，VERA 2B算例布置如圖3所示。具體幾何以及材料信息參見基準(zhǔn)題手冊(cè)[14]。

圖3 VERA 2B俯視圖Fig.3 Top view of VERA 2B

將蒙特卡羅模擬的計(jì)算條件設(shè)置如下：總代數(shù)為500代，其中非活躍代設(shè)置為200代，活躍代設(shè)置為300代，采用ENDF-B/Ⅶ連續(xù)能量點(diǎn)截面數(shù)據(jù)庫(kù)，塊的大小取5，每代模擬的粒子數(shù)為1×105。

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

根據(jù)上述算例設(shè)置，利用NECP-MCX計(jì)算得到VERA 2B算例中keff對(duì)核數(shù)據(jù)的積分靈敏度系數(shù)及69群勒平均靈敏度系數(shù)，并將計(jì)算結(jié)果與MCNP6進(jìn)行對(duì)比。所得keff對(duì)核數(shù)據(jù)的積分靈敏度系數(shù)列于表2。

表2 VERA 2B算例中keff對(duì)核數(shù)據(jù)的積分靈敏度系數(shù)Table 2 Integrated sensitivity coefficients of keff to nuclear data in VERA 2B problem

由表2可發(fā)現(xiàn)，除1H-σelas和1H-σt外，對(duì)于其他核數(shù)據(jù)的靈敏度系數(shù)，NECP-MCX和MCNP6的計(jì)算結(jié)果相對(duì)偏差均小于1%，而1H-σelas和1H-σt的靈敏度系數(shù)計(jì)算結(jié)果相對(duì)偏差分別為-2.99%和-4.12%。這是由于彈性截面的靈敏度系數(shù)計(jì)算涉及到式(13)分子中的第2項(xiàng)和第3項(xiàng)，當(dāng)兩項(xiàng)的量級(jí)相當(dāng)時(shí)，兩項(xiàng)相減會(huì)導(dǎo)致統(tǒng)計(jì)漲落(相對(duì)標(biāo)準(zhǔn)差)較大[15]。MCNP6計(jì)算得到的1H-σelas和1H-σt的靈敏度系數(shù)的相對(duì)標(biāo)準(zhǔn)差也較大，分別為2.90%和4.04%，而NECP-MCX和MCNP6對(duì)這兩個(gè)核素的靈敏度系數(shù)計(jì)算結(jié)果相對(duì)偏差小于MCNP6計(jì)算結(jié)果的3σ，因此可認(rèn)為兩者積分靈敏度系數(shù)計(jì)算結(jié)果符合得較好。

由表2還可發(fā)現(xiàn)，MCNP6計(jì)數(shù)的相對(duì)標(biāo)準(zhǔn)差是NECP-MCX的2.18～2.63倍。這是由于MCNP6使用的是等待法，在當(dāng)前計(jì)算條件設(shè)置下，靈敏度系數(shù)的計(jì)數(shù)次數(shù)為60，而NECP-MCX使用的是重疊塊法，在相同計(jì)算條件設(shè)置下，靈敏度系數(shù)的計(jì)數(shù)次數(shù)為296。根據(jù)式(18)可得到，理論上對(duì)于同一靈敏度系數(shù)的計(jì)算，MCNP6的統(tǒng)計(jì)漲落是NECP-MCX的(296/60)1/2=2.22倍，這一結(jié)果與實(shí)際統(tǒng)計(jì)結(jié)果相符。采用重疊塊法，統(tǒng)計(jì)漲落降低，相應(yīng)的計(jì)算時(shí)間增加，本文未統(tǒng)計(jì)計(jì)算時(shí)間，但根據(jù)文獻(xiàn)[13]結(jié)果，重疊塊法的品質(zhì)因子(FOM)大于等待法。

對(duì)靈敏度系數(shù)分69群統(tǒng)計(jì)，可得到69群勒平均靈敏度系數(shù)，如圖4所示。

對(duì)比圖4中69群勒平均靈敏度系數(shù)可發(fā)現(xiàn)，NECP-MCX的計(jì)算結(jié)果基本上都落在MCNP6計(jì)算結(jié)果的3σ內(nèi)，證明NECP-MCX的靈敏度系數(shù)計(jì)算結(jié)果與MCNP6的吻合較好。同時(shí)可發(fā)現(xiàn)，對(duì)于1H-σelas的勒平均靈敏度系數(shù)結(jié)果，存在較大統(tǒng)計(jì)漲落，同樣是由于在計(jì)算靈敏度系數(shù)時(shí)兩項(xiàng)相減所致，且在能群劃分更細(xì)時(shí)，每個(gè)能群內(nèi)的計(jì)數(shù)更少，靈敏度系數(shù)結(jié)果更難收斂。

圖4 VERA 2B算例中keff對(duì)不同核數(shù)據(jù)的69群勒平均靈敏度系數(shù)Fig.4 69g sensitivity coefficients per unit lethargy of keff to nuclear data in VERA 2B problem

通過(guò)上述對(duì)VERA 2B算例積分靈敏度系數(shù)和69群勒平均靈敏度系數(shù)的計(jì)算結(jié)果對(duì)比，可證明，NECP-MCX計(jì)算得到的keff對(duì)核數(shù)據(jù)的靈敏度系數(shù)結(jié)果與MCNP6吻合較好。

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

本文用NECP-MCX對(duì)AP1000反應(yīng)堆首循環(huán)進(jìn)行建模計(jì)算，反應(yīng)堆處于控制棒全提狀態(tài)，利用NECP-MCX的畫圖功能得到AP1000的堆芯布置圖，如圖5所示。

圖5 AP1000堆芯俯視圖(a)和左視圖(b)Fig.5 Top view (a) and left view (b) of AP1000 core

蒙特卡羅模擬的計(jì)算條件如下：每代粒子數(shù)設(shè)置為2×106，共1 000代，其中前500代為非活躍代，采用ENDF-B/Ⅶ連續(xù)能量點(diǎn)截面數(shù)據(jù)庫(kù)，塊的大小取5。所得AP1000堆芯keff對(duì)連續(xù)能量核數(shù)據(jù)的積分靈敏度系數(shù)列于表3。

由表3可發(fā)現(xiàn)，對(duì)AP1000堆芯keff最敏感的核數(shù)據(jù)分別是235U-ν、235U-σf、235U-σt、238U-σγ、1H-σelas和10B-σt。這是由于235U和238U是首循環(huán)燃料中的主要組成核素，其平均裂變中子數(shù)和裂變截面直接影響中子的產(chǎn)生，對(duì)keff產(chǎn)生直接的影響，同時(shí)AP1000燃料富集度為3%～5%，燃料中238U的比例高達(dá)95%，其俘獲截面較大，會(huì)吸收較多的中子，從而影響keff，且keff的靈敏度系數(shù)為負(fù)。1H是慢化劑H2O中的核素，通過(guò)彈性散射達(dá)到慢化的效果，1H的彈性散射截面越大，其慢化中子的能力越強(qiáng)，反應(yīng)堆中的慢中子越多，發(fā)生裂變反應(yīng)的中子越多，使得keff越大。10B作為硼酸中的主要核素，在反應(yīng)堆中用于控制反應(yīng)性，具有較大的中子吸收截面，因此其總截面靈敏度系數(shù)為負(fù)。

選取表3中積分靈敏度系數(shù)絕對(duì)值較大的6個(gè)反應(yīng)道，得到其69群勒平均靈敏度系數(shù)，如圖6所示。

表3 AP1000算例中keff對(duì)核數(shù)據(jù)的積分靈敏度系數(shù)Table 3 Integrated sensitivity coefficients of keff to nuclear data in AP1000 problem

圖6 AP1000算例中keff對(duì)不同核數(shù)據(jù)的69群勒平均靈敏度系數(shù)Fig.6 69 groups sensitivity coefficients per unit lethargy of keff to nuclear data in AP1000 problem

從圖6可看出，NECP-MCX的計(jì)算結(jié)果與MCNP6符合得較好。靈敏度系數(shù)的形狀與VERA 2B 算例中69群勒平均靈敏度系數(shù)的較近。根據(jù)得到的69群靈敏度系數(shù)以及數(shù)據(jù)庫(kù)中的協(xié)方差矩陣信息，基于“三明治”法則，便可計(jì)算得到核數(shù)據(jù)的不確定度對(duì)AP1000堆芯keff計(jì)算帶來(lái)的不確定度大小，用于后續(xù)的安全分析。

4 總結(jié)

本文基于一階微擾理論推導(dǎo)了keff對(duì)核數(shù)據(jù)的靈敏度系數(shù)公式，采用反復(fù)裂變幾率方法在蒙特卡羅的前向模擬中求解共軛通量，無(wú)需執(zhí)行共軛計(jì)算。采用稀疏矩陣的存儲(chǔ)方式和重疊塊法分別達(dá)到降低存儲(chǔ)內(nèi)存和降低統(tǒng)計(jì)漲落的效果。基于上述理論和方法，在NECP-MCX上完成了連續(xù)能量核數(shù)據(jù)敏感性分析功能模塊的開發(fā)，并對(duì)VERA 2B算例進(jìn)行了驗(yàn)證。結(jié)果表明，NECP-MCX的靈敏度系數(shù)計(jì)算結(jié)果與MCNP6的吻合較好。最后基于AP1000堆芯進(jìn)行了keff對(duì)核數(shù)據(jù)的敏感性分析計(jì)算，得到了對(duì)keff最敏感的核數(shù)據(jù)，分別是235U-ν、235U-σf、235U-σt、238U-σγ、1H-σelas和10B-σt。根據(jù)得到的靈敏度系數(shù)和數(shù)據(jù)庫(kù)中已有的協(xié)方差矩陣信息，基于“三明治”法則，便可得到核數(shù)據(jù)的不確定度對(duì)AP1000堆芯keff計(jì)算導(dǎo)致的不確定度，用于后續(xù)的安全分析。