張 良,江新標,陳立新,楊 寧,張信一
(1.西北核技術研究所,陜西 西安 710024;2.強脈沖輻射環境模擬與效應國家重點實驗室,陜西 西安 710024)
中子價值Φ+(r,E,Ω)的定義[1-2]是:在臨界反應堆中,在r處投入一個能量為E、運動方向為Ω 的中子引起的對穩定功率的貢獻,其物理意義是中子對反應堆裂變功率的貢獻大小。中子價值也稱中子重要性或伴隨中子通量,其概念在1945年由Wigner[3]在反應性方程中首次提出,它在反應堆微擾、敏感性分析和堆芯動態計算方面有著重要應用[4]。Bell和Glasstone[5]在1976年給出了伴隨中子通量的玻耳茲曼輸運方程,中子價值可直接通過求解該方程得到,但在一些幾何較復雜的堆芯中,求解較為困難[6],此時,在滿足擴散近似條件的情況下,可通過求解中子價值的擴散方程得到。在一些蒙特卡羅方法計算動態參數的研究中,也提出了中子價值的多種近似或替代方法,例如Nauchi等[7-8]提出了用中子裂變產生的下一代中子數、Meulekamp等[9]提出了用中子下一代的誘發裂變數來近似代替中子價值,美國洛斯阿拉莫斯國家實驗室的MCNP5-1.60 版本采用中子經若干代后產生的裂變中子總數來近似中子價值函數[10],這些近似方法用于計算動態參數時均得到了較好的結果。
本文采用柵元計算程序WIMS 和擴散計算程序CITATION,計算西安脈沖堆的中子價值分布,研究中子價值隨空間和能量的變化規律,并分析緩發中子有效份額與緩發中子份額、中子代時間與瞬發中子壽命存在差異的原因。
首先采用WIMS程序計算西安脈沖堆的6群群常數,然后采用CITATION 程序進行全堆芯擴散計算,得到堆芯的中子價值分布。
西安脈沖堆堆芯計算模型如圖1所示。脈沖堆堆芯是正六邊形模型,本文為了適應CITATION的幾何描述,在左上和右下角加了兩塊黑體,相當于真空,中子進入黑體后不再返回堆芯。穩態控制棒共5根,最右邊兩根是安全棒,運行中安全棒均被提升至頂位,由于安全棒下部是燃料跟隨體(與燃料棒材料相同),在計算中視為燃料棒。

圖1 西安脈沖堆堆芯模型Fig.1 Core model of Xi'an Pulsed Reactor
WIMS程序應用中子輸運理論來進行柵元計算。
中子輸運方程[2]:

式中:φ(r,E,Ω)為中子在相空間點(r,E,Ω)處的中子通量密度;Σt(r,E)為中子在相空間點(r,E)處的裂變截面;Σs(r,E,Ω→E′,Ω′)為從(E,Ω)到(E′,Ω′)的散射截面;Sf(r,E,Ω)為裂變中子產生率;k為反應堆有效增殖因數。
WIMS程序采用SN 方法對各柵元求解式(1),得到φ(r,E),然后進行群常數計算,本文采用的6群分群結構及說明列于表1。
第g 群的群常數Σm,g可表示為:

式中,m 代表群常數的類型,包括總截面Σt、裂變截面Σf、群轉移截面Σg-g′等。

表1 6群分群結構Table 1 Six-group structure
CITATION 是采用有限差分擴散理論的堆芯分析程序,在堆芯計算中有著廣泛的應用。該程序可處理一維、二維、三維擴散計算問題以及燃耗計算問題,計算精度高,且適用于多種幾何形狀的堆芯。
在由WIMS計算出群常數后,采用CITATION 求解多群共軛擴散方程[2]:

式中:Dg為群擴散系數;ν 為每次裂變產生的中子數。通過求解該方程可得到中子價值(r,E)的分布。本文的計算采用的是二維幾何模型,如圖1 所示,認為堆芯在高度Z 方向是均勻的。
CITATION 程序可同時計算中子通量和中子價值,通過將其計算出的通量分布和連續能量截面的蒙特卡羅程序MCNP[11]計算出的通量分布進行比較,來驗證CITATION 程序用于西安脈沖堆擴散計算的正確性。如圖1所示,選取從左往右第14列,第7~21行的柵元(自上而下編號從1到15,1、14、15號柵元為石墨,2~4、6、10~13號柵元為燃料棒、5號柵元為控制棒,7~9號柵元為水),分別比較第2群(快群)和第6群(熱群)在這15個柵元中的通量分布,結果如圖2所示。兩種程序的計算結果基本符合一致,說明CITATION 程序可較好地用于西安脈沖堆的擴散計算。

圖2 MCNP和CITATION 計算得到的中子通量密度分布對比Fig.2 Comparison of neutron flux density calculated by MCNP and CITATION
如圖1所示,同樣選取從左往右第14列,第7~21行的柵元,分別比較第2群(快群)和第6群(熱群)的中子在這15個柵元中的價值分布,結果如圖3所示。

圖3 中子價值的空間分布Fig.3 Space distribution of neutron importance
對熱群中子,中子價值在2、3號燃料柵元中先增長,然后在4號燃料柵元中有所下降,原因是4號柵元更靠近5號控制棒。然后經歷在控制棒柵元中的急劇降低和在6號燃料柵元中的回升后,在7、8、9號水柵元中先降低再回升,原因是8號水柵元最遠離燃料棒,其中的熱中子最難引起裂變反應,因此中子價值較低。之后中子價值在11號燃料柵元中達到最大,然后依次降低,這是由于12、13號燃料棒逐漸遠離堆芯。快中子的價值整體變化平緩,因為快中子的核反應截面較小,不會因劇烈吸收或引起裂變導致價值突然降低或升高。
分別選擇1、4、5、8號柵元,分別為石墨、燃料、控制棒和水柵元,研究中子價值在這4種柵元中的中子價值隨能量的變化。由圖4可知,石墨柵元的中子價值整體相對較低,原因是石墨柵元主要布置在堆芯外圍。燃料柵元的中子價值隨能量的降低而增大,這是因為燃料的裂變截面隨中子能量的降低而升高。控制棒柵元的中子價值在能量降低時迅速下降,原因是控制棒材料對低能中子的吸收截面很大。水柵元的中子價值隨能量的降低緩慢降低,原因是水對低能中子的吸收較強,而快中子能量更高,有更大的可能性經慢化后進入燃料柵元引起裂變。

圖4 各種柵元在各能群中的中子價值Fig.4 Neutron importance via energy group in different cells
緩發中子有效份額βeff和中子代時間Λ 是反應堆動態特性分析中常用到的重要動態參數,它們的定義式中均考慮了中子價值的影響,βeff和Λ 對應的不考慮中子價值影響的參數分別為緩發中子份額βo 和瞬發中子壽命lP,本文分析βeff與βo、Λ 與lP存在差異的原因。
緩發中子有效份額βeff的定義式為:

中子代時間Λ 的定義式為:

中子代時間Λ 的物理意義是:由φ+(r,E,Ω)作權重的中子數和由φ+(r,E,Ω)作權重的裂變中子產生率之比。即由φ+(r,E,Ω)作權重、每產生一代中子所需的時間。
瞬發中子壽命lP的定義式則為:

對于西安脈沖堆,lP≈93μs,Λ≈36μs。Λ和lP存在明顯的差異,引起這種差異的最重要原因是燃料柵元中的中子價值較其他柵元中的大。如式(5)所示,Λ 定義式中分子中的中子價值包括了所有柵元的中子價值,而分母中的中子價值則只是燃料柵元中的中子價值。而lP的定義式中未考慮中子價值的影響,因此有Λ<lP。
Λ 和lP的相對大小,主要是由燃料柵元的平均中子價值和全堆芯的平均中子價值的相對大小決定的,亦即主要是由中子價值的空間分布決定的。如對裸堆Godiva和Jezebel[12],燃料的中子價值實際上就是全堆芯的中子價值,因此二者的Λ 和lP很接近,相對差別僅為5%左右[13],而這個微小的差別是由于中子價值隨能量的變化引起的,這也同時說明,只考慮中子價值隨能量分布的因素,并不會引起Λ 和lP有很大的差異。
事實上,由于大多類型的反應堆中,燃料柵元的中子價值均大于或接近全堆芯的平均中子價值,因此,一般有Λ≤lP。
本文采用柵元程序WIMS 計算了西安脈沖堆6群群常數,采用擴散程序CITATION 進行全堆芯擴散計算,得到了西安脈沖堆的中子價值在空間和能量上的分布,并獲得以下分析研究結果:
1)得到了中子價值隨空間的變化規律。在選取的15個柵元中,燃料柵元的中子價值較高,控制棒和水柵元中的中子價值相對較低。燃料柵元在遠離堆芯或接近控制棒時,中子價值會降低。快中子價值在各柵元中的變化較熱中子的價值變化平緩。
2)得到了各柵元的中子價值隨能量的變化規律。燃料柵元的中子價值隨能量的降低而增大,原因是燃料發生裂變的可能性隨中子能量的降低而升高。控制棒柵元中的中子價值隨能量的降低而降低,原因是控制棒材料對低能中子有很強的吸收作用。水柵元的中子價值同樣隨能量的降低有緩慢降低,原因是水對低能中子的吸收較強,高能中子則有更大的可能性經慢化后引發裂變。
3)分析了緩發中子有效份額βeff和緩發中子份額βo 存在差異的原因。βeff和βo 的相對大小,取決于中子價值隨能量的分布。對于西安脈沖堆,緩發中子較瞬發中子能量低,其平均價值相對大,因此有βeff>βo。
4)分析了中子代時間Λ 和瞬發中子壽命lP存在差異的原因。Λ 和lP的相對大小,主要是由中子價值隨空間的分布決定的,對于包含西安脈沖堆在內的大多反應堆,燃料柵元中的中子價值較其他柵元中的高,因此有Λ≤lP。
[1] 杜書華.輸運問題的數值模擬[M].湖南:湖南科技出版社,1988.
[2] 謝仲生,鄧力.中子輸運理論數值計算方法[M].西安:西北工業大學出版社,2005.
[3] WIGNER E P.The collected works of Eugene Paul Wigner:Effect of small perturbations on pile period[M].Berlin:Springer Verlag,1992.
[4] HENRY A F.Nuclear reactor analysis[M].Massachusetts:MIT Press,1975.
[5] BELL G I,GLASSTONE S.Nuclear reactor theory[M].New York:Van Nostrand Reinhold Company,1970.
[6] FEGHHI S A H,SHAHRIARI M,AFARIDEH H.Calculation of neutron importance function in fissionable assemblies using Monte Carlo method[J].Annals of Nuclear Energy,2007,34:514-520.
[7] NAUCHI Y,KAMEYAMA T.Proposal of direct calculation of kinetic parametersβeffand Λ based on continuous energy Monte Carlo method[J].Journal of Nuclear Science and Technology,2005,42(6):503-514.
[8] NAUCHI Y,KAMEYAMA T.Development of calculation technique for iterated fission probability and reactor kinetic parameters using continuous-energy Monte Carlo method[J].Journal of Nuclear Science and Technology,2010,47(11):977-990.
[9] MEULEKAMP R K,van der MARCK S C.Calculating the effective delayed neutron fraction using Monte Carlo[J].Nuclear Science and Engineering,2006,152:142-148.
[10]KIEDROWSKI B C,BOOTH T E,BROWN F B,et al.MCNP5-1.60feature enhancements &manual clarifications[R].USA:Los Alamos National Laboratory,2010.
[11]BRIESMEISTER J F.Mont Carlo N-particle transport[R].USA:Los Alamos National Laboratory,2000.
[12]OECD/NEA Nuclear Science Committee.International handbook of evaluated criticality safety benchmark experiments[M].USA:OECD/NEA,2008.
[13]KIEDROWSKI B C,BROWN F.Adjoint wighted kinetics parameters with continuous energy Monte Carlo[J].Tran Am Nucl Soc,2009,100:297-299.