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

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

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

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

1 理論方法

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

1.1 一階微擾理論

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

(1)

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

輸運-碰撞項算子B表示為：

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

(2)

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

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

(3)

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

(4)

共軛輸運-碰撞項算子B*表示為：

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

(5)

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

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

(6)

引入擾動：

F′=F+δF

B′=B+δB

k′=k+δk

ψ′=ψ+δψ

(7)

則擾動后的中子輸運方程表示為：

(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)

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

(14)

散射項：

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

(15)

碰撞項：

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

(16)

分母表示為：

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

(17)

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

1.2 反復裂變幾率

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

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

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

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

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

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

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

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

由于蒙特卡羅模擬中隨機變量的統(tǒng)計標準差σ滿足：

(18)

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

2 程序驗證

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

2.1 算例設(shè)置

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

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

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

2.2 計算結(jié)果

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

表2 VERA 2B算例中keff對核數(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外，對于其他核數(shù)據(jù)的靈敏度系數(shù)，NECP-MCX和MCNP6的計算結(jié)果相對偏差均小于1%，而1H-σelas和1H-σt的靈敏度系數(shù)計算結(jié)果相對偏差分別為-2.99%和-4.12%。這是由于彈性截面的靈敏度系數(shù)計算涉及到式(13)分子中的第2項和第3項，當兩項的量級相當時，兩項相減會導致統(tǒng)計漲落(相對標準差)較大[15]。MCNP6計算得到的1H-σelas和1H-σt的靈敏度系數(shù)的相對標準差也較大，分別為2.90%和4.04%，而NECP-MCX和MCNP6對這兩個核素的靈敏度系數(shù)計算結(jié)果相對偏差小于MCNP6計算結(jié)果的3σ，因此可認為兩者積分靈敏度系數(shù)計算結(jié)果符合得較好。

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

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

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

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

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

3 AP1000計算結(jié)果

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

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

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

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

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

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

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

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

4 總結(jié)

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