連續能量核數據敏感性分析及其在AP1000中的應用

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

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

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

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

1 理論方法

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

1.1 一階微擾理論

穩態中子輸運方程可表示為：

(1)

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

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

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

(2)

產生項算子F表示為：

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

(3)

同樣，穩態共軛中子輸運方程可表示為：

(4)

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

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

(5)

共軛產生項算子F*表示為：

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

(6)

引入擾動：

F′=F+δF

B′=B+δB

k′=k+δk

ψ′=ψ+δψ

(7)

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

(8)

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

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

(9)

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

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

(10)

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

(11)

靈敏度系數的定義為：

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

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

1.2 反復裂變幾率

從式(14)～(17)可看出，求解連續能量核數據的靈敏度系數需得到共軛通量以及核反應率，統計核反應率在蒙特卡羅程序中易實現，因此求解靈敏度系數的關鍵在于求解共軛通量。

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

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

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

表1 核反應率計數數組非零值比例Table 1 None-zero value rate in reaction rate array

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

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

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

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

(18)

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

2 程序驗證

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

2.1 算例設置

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

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

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

2.2 計算結果

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

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

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

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

對靈敏度系數分69群統計，可得到69群勒平均靈敏度系數，如圖4所示。

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

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

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

3 AP1000計算結果

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

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

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

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

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

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

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

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

4 總結

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