一種采用蒙卡程序計算兩群中子價值分布的方法①

安偉健 趙澤昊 霍紅磊

(中國原子能科學研究院反應堆工程技術研究部 北京 102413)

在開展反應堆物理特性分析時，中子價值往往可提供重要的幫助。中子價值的求解在確定論程序中較為容易，只需對中子輸運方程的共軛方程進行求解，其求解過程和中子通量密度的求解同樣容易。然而，對于蒙特卡羅程序（簡稱蒙卡程序），中子價值的求解則相對較為困難。

中子價值的定義為：臨界系統(tǒng)中，在某一位置投入某一能量和運動方向的中子，由該中子引起的對系統(tǒng)穩(wěn)定功率的貢獻，即為該中子的價值。美國的Hurwitz還指出：臨界系統(tǒng)中，由1個中子引起的每代裂變次數(shù)在代數(shù)趨于無窮時（此時中子分布已趨于穩(wěn)定）會趨于一定值，即為中子價值。根據(jù)該理論，國際上針對蒙卡程序陸續(xù)開發(fā)了幾種中子價值計算方法，分別將反復裂變幾率、下一代裂變幾率或下一代裂變中子數(shù)目的統(tǒng)計結果作為中子價值；在國內，汪量子等曾以反復裂變幾率作為中子價值，將該方法融入到MCNP以及多群蒙卡程序MCMG中，獲得了很好的計算效果[1]。

本文基于對中子價值物理意義的認識，提出了一種采用蒙卡程序計算兩群中子價值分布的新方法，并將該方法成功應用于實際的反應堆物理特性分析工作中，取得了較好的效果。

1 方法描述

文獻[2]提出了一種求解燃料內兩群中子平均價值的比值的方法。該方法基于對反應堆內中子行為的認識，畫出了燃料內熱群中子的輸運流程圖，見圖1。圖中，P0為燃料內熱群中子的首次飛行逃脫幾率，C為燃料內熱群中子在碰撞時發(fā)生散射的幾率，v2為熱群中子的平均裂變中子數(shù)；Σ2f為熱群中子的宏觀裂變截面； Σ2a為熱群中子的宏觀吸收截面。

根據(jù)中子價值守恒原理，初始中子的價值應當?shù)扔谄浜蟠凶拥目們r值。圖中，初始中子的后代中子可分為兩部分：在燃料里被吸收產生新一代的中子，以及泄漏至燃料外的中子。因而圖1可簡化為圖2的形式，并可建立關系式（1）。其中，PFF,2為燃料內熱群中子在燃料內被吸收的幾率，PFO,2為燃料內熱群中子泄漏至燃料外的幾率；和分別為燃料內熱群及快群中子的平均價值，為燃料外熱群中子的平均價值。由此即可得到燃料內外兩區(qū)的兩群中子平均價值之間的關系（對于非臨界系統(tǒng)，keff不能省略，用以維持不同代的中子價值守恒）。

上述分析方法可進行拓展：若將反應堆劃分為多個區(qū)，對各相鄰區(qū)域分別開展分析，即可獲得各區(qū)域的兩群中子平均價值之間的關系式，進而可求得整個反應堆內的兩群中子價值分布情況。以圖3模型為例進行說明，將反應堆由內至外劃分為三區(qū)，其中活性區(qū)分為兩區(qū)（a區(qū)和b區(qū)），反射層為一區(qū)（c區(qū)）。

對a區(qū)的熱群中子進行分析，可畫出其輸運簡圖（見圖4），并建立關系式（2）。其中，va,2為a區(qū)熱群中子的平均裂變中子數(shù)；Σa,f,2為a區(qū)熱群中子的宏觀裂變截面；Σa,a,2為a區(qū)熱群中子的宏觀吸收截面，其余各參數(shù)的命名規(guī)則可參考關系式（1）。

此外，對a區(qū)的快群中子進行分析，也可畫出其輸運簡圖（見圖5），并建立關系式（3）。其中，Pa,12為a區(qū)快群中子慢化至熱群的幾率，其他參數(shù)命名規(guī)則可參考關系式（1）和（2）。

類似的，對b區(qū)的熱群中子進行分析，也可畫出其輸運簡圖（見圖6），并建立關系式（4）。參數(shù)命名規(guī)則可參考關系式（1）和（2）。

同樣，對b區(qū)的快群中子、c區(qū)的熱群及快群中子，均可分別建立對應的關系式如下：

以上各式中的各項系數(shù)均可較容易由蒙卡程序MCNP計算并處理得到。將關系式（2）-（7）聯(lián)立，可得6個關系式和6個未知量，實際計算過程中發(fā)現(xiàn)，該方程組并無法求出各未知量的確定解，因為計算過程中會出現(xiàn)兩個等同的關系式，但可由其中任意5個關系式求出各未知量的相對值，也就是各區(qū)的兩群中子平均價值之間的相對關系。以此為基礎，若將反應堆較細致地劃分為多個區(qū)域，采用該方法就可以獲得整個反應堆內的兩群中子價值分布情況。

圖1 燃料內熱群中子的輸運流程圖

圖2 燃料內熱群中子的輸運簡圖

2 方法應用

對一個采用CERMET燃料（燃料成分為W-60%UO2-6%Gd2O3，包殼為W-25%Re）的空間核反應堆的臨界安全特性進行計算時發(fā)現(xiàn)：在反應堆進水的情況下，隨著進水量的增加，系統(tǒng)的effk呈現(xiàn)出先增大后減小的特性。下面采用上述方法，從中子價值的角度對該現(xiàn)象進行分析。

首先，將該反應堆簡化為球形模型，見圖7，其中活性區(qū)半徑為18cm，BeO反射層厚度為10cm?；钚詤^(qū)進水時（假設水在活性區(qū)內均勻分布），系統(tǒng)effk隨進水量的變化曲線見圖8，在堆內水密度小于約0.03g/cm3時，k eff隨進水量的增加而增大；而在堆內水密度大于約0.03g/cm3時，k eff將隨進水量的增加而減小。

對該現(xiàn)象進行分析：活性區(qū)進水時，對系統(tǒng)的直接影響是增加了中子慢化性能，軟化了堆芯的中子能譜。而中子能量的變化將影響其價值，也就是對effk的貢獻，因而，考慮從中子價值的角度來進行研究，求解出堆內不同能群中子的價值分布，進而評估中子能譜的改變對effk的影響。

圖3 反應堆分區(qū)模型

圖4 a區(qū)熱群中子的輸運簡圖

圖5 a區(qū)快群中子的輸運簡圖

圖6 b區(qū)熱群中子的輸運簡圖

由于該反應堆的中子能譜非常硬（無慢化劑，且燃料內含較多W、Re、Gd熱中子毒物），熱群中子所占份額極小，因此考慮以5.5 keV作為分界能，第二群中子包括共振中子和熱中子，第一群中子則包括快中子和中能中子，裂變中子均為第一群中子。

將該球形反應堆模型按同心圓分成23區(qū)，其中活性區(qū)16區(qū)，反射層7區(qū)，從里到外用字母進行編號，見圖9。在靠近活性區(qū)與反射層的交界處劃分較細，這是由于該處中子通量密度隨空間變化較快。

采用本文所提方法，對每個區(qū)的兩群中子平均價值均可建立一個線性方程，這樣一共可建立46個線性方程。而每個區(qū)均有兩個未知量（也就是兩群中子的平均價值），因此一共有46個未知量。在實際計算過程中會出現(xiàn)兩個等同的關系式，可任選其一。這樣，由45個方程可求解出46個未知量的相對值。令活性區(qū)最外層的第二群中子平均價值為1，將該線性方程組寫為矩陣形式，可用MATLAB程序很容易求解。對于不同進水量的情況，計算結果如圖10～12所示。

圖7 反應堆簡化球形模型

圖8 反應堆 effk 隨進水量的變化

圖9 分區(qū)的計算模型

根據(jù)圖10，活性區(qū)未進水時，在距離球心約0～8cm范圍內，第一群中子的價值大于第二群中子，而在距球心8～18cm的區(qū)域（占活性區(qū)絕大部分區(qū)域），第一群中子的價值小于第二群中子，這使得在活性區(qū)進水導致能譜軟化時，8～18cm的區(qū)域對反應性的正貢獻占主導地位，從而使effk上升。再由圖11和圖12可以看出，隨著進水量的增加，兩群中子價值分布曲線的交叉點有明顯的右移趨勢。到水密度為0.1g/cm3時，活性區(qū)內絕大部分區(qū)域的第一群中子價值均大于第二群中子。這使得此時隨著進水量的增加，keff將呈現(xiàn)下降趨勢。這就解釋了圖8中effk隨進水量的增加呈現(xiàn)出先增大后減小的現(xiàn)象。

圖10 堆芯兩群中子價值相對分布（未進水）

圖11 堆芯兩群中子價值相對分布（水密度0.03 g/cm3）

基于上述計算和分析結果，從物理方面可對該現(xiàn)象可作如下解釋：在活性區(qū)靠近中心的區(qū)域，對中子價值起主導作用的是中子利用系數(shù)，由于W、Re、Gd等第二群中子毒物的存在，使得第二群中子的利用系數(shù)低于第一群中子，因此其價值低于第一群中子；在活性區(qū)靠外的區(qū)域，泄漏是影響中子價值的主導因素，第一群中子的泄漏率要比第二群中子高，導致第二群中子的價值高于第一群中子。在堆內未進水或進水量較少時，隨著進水量的增加，中子的利用系數(shù)減小，不泄漏幾率增大，后者對系統(tǒng)的影響大于前者，使得effk增大；在堆內進水量較多時，隨著進水量的增加，中子的利用系數(shù)減小，不泄漏幾率增大，此時前者占主導地位，使得effk減小。

3 結果驗證

采用微擾方法對上述中子價值計算結果進行驗證。具體如下：

CERMET燃料里包含一定量的Gd2O3。若對Gd的量作微小的調整，對中子通量密度的分布影響很小，可采用微擾方法計算反應性的變化量。該值也可通過計算調整前后的effk變化量來求得。若上述中子價值的計算結果正確，則這兩種方法的計算結果應當相符。

圖12 堆芯兩群中子價值相對分布（水密度0.1g/cm3）

將CERMET燃料中的Gd減少10%。按照微擾理論[3]，反應性的變化量應當?shù)扔谶@10%的Gd所吸收的中子價值與系統(tǒng)內中子總價值的比值：

式中，ΣGd,1、ΣGd,2分別為Gd的宏觀吸收截面，RGd,i,1、分別為第i區(qū)內Gd對第一群和第二群中子的吸收率，Ni為第i區(qū)燃料的裂變中子產生率。各參數(shù)均可由MCNP程序很容易算得。結合之前的計算結果，可算得由于Gd的減少所引起的反應性變化量為：另一方面，通過計算調整前后的effk變化量也可得到Δρ的數(shù)值，該方法的計算結果為：兩種方法的計算結果差別僅為4.2%。這說明該方法的計算結果是合理可信的。

4 結語

針對一般蒙卡程序無法計算中子價值的問題，提出了一種采用蒙卡程序計算兩群中子價值分布的新方法，并將該方法成功應用于實際反應堆的物理特性分析中。所存在的不足之處在于：由于手頭上沒有其他有效的中子價值計算程序，因此沒能對中子價值分布的計算結果進行直接校驗，而是通過微擾方法進行了間接的驗證。