朱金輝 黃流興 牛勝利

（西北核技術(shù)研究所 西安 710024）

蒙特卡羅方法可用于超臨界系統(tǒng)的裂變增值模擬，大型蒙特卡羅模擬程序MCNP在反應(yīng)堆臨界分析計算中應(yīng)用廣泛[1?5]。對于超臨界系統(tǒng)，由于其中子不斷增值，在進行蒙特卡羅中子輸運模擬時通常需用特殊的處理方法。在MCNP中，采用裂變驛站方法[6]，通過增值系數(shù)作為校正因子以控制模擬中子數(shù)的增加，可計算給出系統(tǒng)的增值系數(shù)和中子平均去除時間等參數(shù)。在超臨界系統(tǒng)分析中，反映系統(tǒng)參數(shù)隨時間變化的增值時間常數(shù)，是超臨界系統(tǒng)分析中的關(guān)鍵參數(shù)，增值時間常數(shù)與增值系數(shù)的關(guān)系也受到重點關(guān)注[7]。為反映系統(tǒng)增值的時間特性而開發(fā)的能模擬中子g耦合輸運問題的粒子輸運蒙卡模擬程序TOPAN，其中子輸運模塊根據(jù)ENDF/B中子截面庫開發(fā)，具備裂變驛站、時間驛站、粒子標識等功能，可模擬超臨界系統(tǒng)中的中子輸運，給出超臨界系統(tǒng)的增值系數(shù)、增值時間常數(shù)等。

本文用TOPAN程序?qū)ΤR界系統(tǒng)進行模擬計算，給出系統(tǒng)的增值系數(shù)和中子平均去除時間，并與MCNP結(jié)果進行對比，驗證其超臨界系統(tǒng)計算準確性和模擬鏈式反應(yīng)過程的可靠性。用TOPAN程序的時間驛站功能計算超臨界系統(tǒng)的增值時間常數(shù)，并分析其與增值系數(shù)和平均裂變時間的關(guān)系。

1 系統(tǒng)增值系數(shù)計算

TOPAN程序中，用裂變驛站法可由增值系數(shù)的直接定義算得系統(tǒng)的增值系數(shù)keff，即

其中Ni是第i代裂變產(chǎn)生的中子總數(shù)。在每一代裂變中子輸運過程中，裂變中子不再繼續(xù)輸運跟蹤，而是存儲為驛站源。在下一裂變代中子輸運前，將裂變驛站源抽樣得到中子輸運的源粒子。經(jīng)最初n代裂變中子輸運后，形成較為穩(wěn)定的中子分布參數(shù)。在后續(xù)裂變代的輸運中根據(jù)每一裂變代源中子的總權(quán)重和產(chǎn)生的裂變中子的總權(quán)重按式(1)給出系統(tǒng)增值系數(shù)keff。

本文用一個钚球(半徑r1=5 cm)外包一個鎢球殼反射層(半徑r2分別取5、10、15和20 cm)的簡單模型。用TOPAN程序和MCNP程序計算其特征值增值系數(shù)keff，結(jié)果見表1。計算給出的增值系數(shù)keff與MCNP程序基本一致，證明TOPAN程序在鏈式增值反應(yīng)中子輸運模擬計算的可靠性。而有反射層的模型的增值系數(shù)明顯大于無反射層的模型，且增值系數(shù)隨著反射層厚度增加，但增值系數(shù)的增長速率隨反射層厚度的增加而變慢。

2 系統(tǒng)中子平均去除時間計算

中子平均去除時間(逃逸、吸收或裂變)是超臨界系統(tǒng)分析中的重要參數(shù)。TOPAN程序中用式(2)統(tǒng)計平均去除時間。

其中We和Te為逃逸粒子的權(quán)重和時間，Wc和Tc為被吸收粒子的權(quán)重和時間，Wf和Tf為發(fā)生裂變反應(yīng)粒子的權(quán)重和時間，采用與§1相同的計算模型，計算給出的平均去除時間tr如表2所示。計算

表2 平均去除時間tr的計算結(jié)果*Table 2 The results of the prompt removal lifetime

得到的平均去除時間與MCNP程序計算結(jié)果相當(dāng)，偏差小于 2%。同時的增加，平均去除時間隨反射層厚度迅速增加，可見反射層的存在使中子輸運和中子逃逸時間大大增加。

系統(tǒng)中子去除有三種途徑：逃逸、裂變和吸收，其中裂變發(fā)生的時間對系統(tǒng)增值的時間特性有直接影響。TOPAN程序可以求出三種不同去除途徑的平均時間及比例。定義平均逃逸時間平均逃逸時間te，平均裂變時間tf，平均吸收時間tc如下：

計算結(jié)果見表3。三種機制中，平均裂變時間tf最短。隨反射層厚度增加，平均逃逸時間和平均吸收時間迅速增加，逃逸的權(quán)重比則不斷減小。

表3 不同去除機制的平均時間計算結(jié)果(單位：ns)Table 3 The results of removal lifetime(in ns) for various reactions

3 系統(tǒng)增值時間常數(shù)計算

對于一個超臨界系統(tǒng)，系統(tǒng)增值時間常數(shù)反映其系統(tǒng)參數(shù)（如功率、中子數(shù)）隨的時間變化。利用程序的時間驛站功能，可給出系統(tǒng)的增值時間常數(shù)l?；谇笆瞿Ｐ停O(shè)14.06 MeV中子源位于系統(tǒng)中心，源強為8′105，取時間步長為 0.1 ns，粒子輸運到時間步結(jié)束時間時存儲為時間驛站源，在下一時間開始，從時間驛站源中抽樣得到一定數(shù)量的源粒子開始下一時間的輸運計算。由于r1=r2=5 cm的系統(tǒng)接近臨界，將其改為r1=r2=6 cm的模型，同時取r1=5 cm，r2分別為10、15和20 cm，計算得到系統(tǒng)的功率隨時間變化的關(guān)系如圖1(a)，系統(tǒng)內(nèi)中子數(shù)隨時間的變化關(guān)系如圖1(b)。

圖1 系統(tǒng)功率(a)和中子總數(shù)(b)隨時間的變化Fig. 1 The system power (a) and neutron counts (b) as function of the time.

由圖1，經(jīng)過最初的時間0t(10–20 ns)后，系統(tǒng)的功率釋放和中子總數(shù)呈指數(shù)增長，即：

系統(tǒng)中的中子數(shù)隨時間的變化可表示為：

實際上，系統(tǒng)功率增值和中子增值的時間常數(shù)l值基本一致，計算給出的結(jié)果列于表4。

表4 超臨界系統(tǒng)增值相關(guān)參數(shù)計算結(jié)果Table 4 The multiplication parameters of the supercritical system.

4 超臨界系統(tǒng)增值時間常數(shù)與增值系數(shù)的關(guān)系

顯然，超臨界系統(tǒng)的增值時間常數(shù)l與增值系數(shù)effk直接相關(guān)。為將上述系統(tǒng)的增值時間常數(shù)與其增值系數(shù)聯(lián)系起來，可先考慮一種理想情況，即一代裂變中子從產(chǎn)生到發(fā)生下一代裂變反應(yīng)的時間分布滿足單一負指數(shù)分布，即

則系統(tǒng)每過 1/時間，系統(tǒng)增值 exp(k?1)倍，則系統(tǒng)平均裂變時間和增值時間常數(shù)值及增值系數(shù)的關(guān)系為

上述分析表明，系統(tǒng)增值時間常數(shù)主要受增值系數(shù)keff和平均裂變時間的影響。然而，式(8)的前提條件—中子裂變時間分布為單一負指數(shù)分布—與實際情況并不相符，裂變發(fā)生的時間分布可近似表示為式(9)的多指數(shù)分布

為此可參照式(8)，根據(jù)系統(tǒng)的增值時間常數(shù)和增值系數(shù)定義等效裂變時間為

計算給出的等效裂變時間列于表4。對于不帶反射層的系統(tǒng)，其等效裂變時間與平均裂變時間相當(dāng)，略小于系統(tǒng)的平均裂變時間；隨著放射層厚度的增加，平均裂變時間與等效裂變時間的差值逐漸增大。原因簡單分析如下：裂變發(fā)生時間分布式(9)中既有下降較快的快成分，也有下降較慢的慢成分；慢成分的裂變發(fā)生時間（正比于較?。┹^長，使平均裂變時間增加，而對系統(tǒng)增值貢獻（正比于）為較??；綜合起來分析，系統(tǒng)的平均裂變時間大于等效裂變時間；對于具有較厚反射層的系統(tǒng)，其慢成分的發(fā)生時間可能長達數(shù)百ns甚至幾個μs，大大增加了系統(tǒng)的平均裂變時間，從而使系統(tǒng)的平均裂變時間和等效裂變時間差值變得較大。

表4數(shù)據(jù)還表明，并不一定增值系數(shù)keff大的系統(tǒng)其系統(tǒng)功率的增值速度就越大，分析計算時還應(yīng)考慮平均裂變時間對系統(tǒng)功率增值的影響。

5 小結(jié)

本文建立的蒙卡程序能夠有效用于超臨界系統(tǒng)中子輸運模擬，能夠給出系統(tǒng)增值系統(tǒng)、增值時間常數(shù)等參數(shù)。對于超臨界系統(tǒng)，其系統(tǒng)增值時間常數(shù)除和增值系數(shù)keff有關(guān)以外，還主要受平均裂變時間的影響。實際超臨界系統(tǒng)中的等效裂變時間一般小于平均裂變時間，隨著反射層厚度增加，兩者差別隨之增大。

本文的計算方法可用于超臨界系統(tǒng)增值反應(yīng)相關(guān)時間特性分析，對反應(yīng)堆啟動和脈沖堆運行等具有一定的理論意義和應(yīng)用價值。

1 李樹, 謝仲生, 鄧力, 等. 蒙特卡羅方法在反應(yīng)堆物理計算中的應(yīng)用[J], 核科學(xué)與工程, 2003, 23(2), 188-192 LI Shu, XIE Zhongsheng, DENG Li,et al. Application of Monte Carlo method in the calculation of nuclear reactor physics[J], Chinese Journal of Nuclear Science and Engineering, 2003, 23(2), 188-192

2 鐘兆鵬, 施工, 胡永明. MCNP程序在反應(yīng)堆臨界計算中的應(yīng)用[J], 核動力工程, 2003, 24(1), 8-11 ZHONG Zhaopeng, SHI Gong, HU Yongming.Application of MCNP in the criticality, calculation for reactors[J], Nuclear Power Engineering, 2003, 24(1),8-11

3 邱小平, 黎學(xué)川, 王建華, 次臨界堆芯參數(shù)變化對Keff值的影響[J], 南華大學(xué)學(xué)報(自然科學(xué)版), 2005, 19(1),42-46 QIU Xiaoping, LI Xuechuan, WANG Jianhua, The Influence onKeffby parameter change in subcritical reactor[J], Journal of Nanhua University(Science and Technology), 2005, 19(1), 42-46

4 邱有恒, 鄧力. 次臨界a本征值計算中的時間選取方法[J], 核科學(xué)與工程, 2007, 27(2): 182-186 QIU Youheng, DENG Li. The method of selection time interval in the calculation of subcritical a eigenvalue[J],Chinese Journal of Nuclear Science and Engineering,2007, 27(2): 182-186

5 王瑞宏, 黃正豐, 江松, 等. 計算a本征值的蒙特卡羅期望估計方法[J], 核科學(xué)與工程, 2008, 28(3): 204-209 WANG Ruihong , HUANG Zhengfeng , JIANG Song,et al. Monte Carlo expected estimation method foraeigenvalue calculation[J], Chinese Journal of Nuclear Science and Engineering, 2008, 28(3): 204-209

6 Judith F. Briesmeister, Editor. MCNP-A General Monte Carlo N-particle Transport Code. Version 4B[M]. Los Alamos National Laboratory, LA-12625-M , 1997

7 邱有恒, 李茂生, 鄧力, 次臨界系統(tǒng)中a本征值與k本征值的關(guān)系[J], 清華大學(xué)學(xué)報(自然科學(xué)版), 2007,47(S1): 967-971 QIU Youheng, LI Maosheng, DENG Li, Relation ship between theaeigenvalue and thekeigenvalue in subcritical systems[J], J Tsinghua Univ (Sci & Tech) ,2007, 47(S1): 967-971