張 良，江新標(biāo)，陳立新，楊 寧，張信一

（1.西北核技術(shù)研究所，陜西 西安 710024；2.強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710024）

中子價(jià)值Φ＋（r，E，Ω）的定義［1-2］是：在臨界反應(yīng)堆中，在r處投入一個(gè)能量為E、運(yùn)動(dòng)方向?yàn)棣?的中子引起的對(duì)穩(wěn)定功率的貢獻(xiàn)，其物理意義是中子對(duì)反應(yīng)堆裂變功率的貢獻(xiàn)大小。中子價(jià)值也稱中子重要性或伴隨中子通量，其概念在1945年由Wigner［3］在反應(yīng)性方程中首次提出，它在反應(yīng)堆微擾、敏感性分析和堆芯動(dòng)態(tài)計(jì)算方面有著重要應(yīng)用［4］。Bell和Glasstone［5］在1976年給出了伴隨中子通量的玻耳茲曼輸運(yùn)方程，中子價(jià)值可直接通過求解該方程得到，但在一些幾何較復(fù)雜的堆芯中，求解較為困難［6］，此時(shí)，在滿足擴(kuò)散近似條件的情況下，可通過求解中子價(jià)值的擴(kuò)散方程得到。在一些蒙特卡羅方法計(jì)算動(dòng)態(tài)參數(shù)的研究中，也提出了中子價(jià)值的多種近似或替代方法，例如Nauchi等［7-8］提出了用中子裂變產(chǎn)生的下一代中子數(shù)、Meulekamp等［9］提出了用中子下一代的誘發(fā)裂變數(shù)來近似代替中子價(jià)值，美國洛斯阿拉莫斯國家實(shí)驗(yàn)室的MCNP5-1.60 版本采用中子經(jīng)若干代后產(chǎn)生的裂變中子總數(shù)來近似中子價(jià)值函數(shù)［10］，這些近似方法用于計(jì)算動(dòng)態(tài)參數(shù)時(shí)均得到了較好的結(jié)果。

本文采用柵元計(jì)算程序WIMS 和擴(kuò)散計(jì)算程序CITATION，計(jì)算西安脈沖堆的中子價(jià)值分布，研究中子價(jià)值隨空間和能量的變化規(guī)律，并分析緩發(fā)中子有效份額與緩發(fā)中子份額、中子代時(shí)間與瞬發(fā)中子壽命存在差異的原因。

1 計(jì)算方法

首先采用WIMS程序計(jì)算西安脈沖堆的6群群常數(shù)，然后采用CITATION 程序進(jìn)行全堆芯擴(kuò)散計(jì)算，得到堆芯的中子價(jià)值分布。

1.1 西安脈沖堆堆芯模型

西安脈沖堆堆芯計(jì)算模型如圖1所示。脈沖堆堆芯是正六邊形模型，本文為了適應(yīng)CITATION的幾何描述，在左上和右下角加了兩塊黑體，相當(dāng)于真空，中子進(jìn)入黑體后不再返回堆芯。穩(wěn)態(tài)控制棒共5根，最右邊兩根是安全棒，運(yùn)行中安全棒均被提升至頂位，由于安全棒下部是燃料跟隨體（與燃料棒材料相同），在計(jì)算中視為燃料棒。

圖1 西安脈沖堆堆芯模型Fig.1 Core model of Xi'an Pulsed Reactor

1.2 WIMS程序計(jì)算6群群常數(shù)

WIMS程序應(yīng)用中子輸運(yùn)理論來進(jìn)行柵元計(jì)算。

中子輸運(yùn)方程［2］：

式中：φ（r，E，Ω）為中子在相空間點(diǎn)（r，E，Ω）處的中子通量密度；Σt（r，E）為中子在相空間點(diǎn)（r，E）處的裂變截面；Σs（r，E，Ω→E′，Ω′）為從（E，Ω）到（E′，Ω′）的散射截面；Sf（r，E，Ω）為裂變中子產(chǎn)生率；k為反應(yīng)堆有效增殖因數(shù)。

WIMS程序采用SN 方法對(duì)各柵元求解式（1），得到φ（r，E），然后進(jìn)行群常數(shù)計(jì)算，本文采用的6群分群結(jié)構(gòu)及說明列于表1。

第g 群的群常數(shù)Σm，g可表示為：

式中，m 代表群常數(shù)的類型，包括總截面Σt、裂變截面Σf、群轉(zhuǎn)移截面Σg-g′等。

表1 6群分群結(jié)構(gòu)Table 1 Six-group structure

1.3 CITATION 程序擴(kuò)散計(jì)算

CITATION 是采用有限差分?jǐn)U散理論的堆芯分析程序，在堆芯計(jì)算中有著廣泛的應(yīng)用。該程序可處理一維、二維、三維擴(kuò)散計(jì)算問題以及燃耗計(jì)算問題，計(jì)算精度高，且適用于多種幾何形狀的堆芯。

在由WIMS計(jì)算出群常數(shù)后，采用CITATION 求解多群共軛擴(kuò)散方程［2］：

式中：Dg為群擴(kuò)散系數(shù)；ν 為每次裂變產(chǎn)生的中子數(shù)。通過求解該方程可得到中子價(jià)值（r，E）的分布。本文的計(jì)算采用的是二維幾何模型，如圖1 所示，認(rèn)為堆芯在高度Z 方向是均勻的。

2 CITATION 擴(kuò)散計(jì)算的正確性驗(yàn)證

CITATION 程序可同時(shí)計(jì)算中子通量和中子價(jià)值，通過將其計(jì)算出的通量分布和連續(xù)能量截面的蒙特卡羅程序MCNP［11］計(jì)算出的通量分布進(jìn)行比較，來驗(yàn)證CITATION 程序用于西安脈沖堆擴(kuò)散計(jì)算的正確性。如圖1所示，選取從左往右第14列，第7～21行的柵元（自上而下編號(hào)從1到15，1、14、15號(hào)柵元為石墨，2～4、6、10～13號(hào)柵元為燃料棒、5號(hào)柵元為控制棒，7～9號(hào)柵元為水），分別比較第2群（快群）和第6群（熱群）在這15個(gè)柵元中的通量分布，結(jié)果如圖2所示。兩種程序的計(jì)算結(jié)果基本符合一致，說明CITATION 程序可較好地用于西安脈沖堆的擴(kuò)散計(jì)算。

圖2 MCNP和CITATION 計(jì)算得到的中子通量密度分布對(duì)比Fig.2 Comparison of neutron flux density calculated by MCNP and CITATION

3 西安脈沖堆中子價(jià)值的計(jì)算結(jié)果及分析

3.1 中子價(jià)值的空間分布

如圖1所示，同樣選取從左往右第14列，第7～21行的柵元，分別比較第2群（快群）和第6群（熱群）的中子在這15個(gè)柵元中的價(jià)值分布，結(jié)果如圖3所示。

圖3 中子價(jià)值的空間分布Fig.3 Space distribution of neutron importance

對(duì)熱群中子，中子價(jià)值在2、3號(hào)燃料柵元中先增長，然后在4號(hào)燃料柵元中有所下降，原因是4號(hào)柵元更靠近5號(hào)控制棒。然后經(jīng)歷在控制棒柵元中的急劇降低和在6號(hào)燃料柵元中的回升后，在7、8、9號(hào)水柵元中先降低再回升，原因是8號(hào)水柵元最遠(yuǎn)離燃料棒，其中的熱中子最難引起裂變反應(yīng)，因此中子價(jià)值較低。之后中子價(jià)值在11號(hào)燃料柵元中達(dá)到最大，然后依次降低，這是由于12、13號(hào)燃料棒逐漸遠(yuǎn)離堆芯?？熘凶拥膬r(jià)值整體變化平緩，因?yàn)榭熘凶拥暮朔磻?yīng)截面較小，不會(huì)因劇烈吸收或引起裂變導(dǎo)致價(jià)值突然降低或升高。

3.2 中子價(jià)值隨能量的分布

分別選擇1、4、5、8號(hào)柵元，分別為石墨、燃料、控制棒和水柵元，研究中子價(jià)值在這4種柵元中的中子價(jià)值隨能量的變化。由圖4可知，石墨柵元的中子價(jià)值整體相對(duì)較低，原因是石墨柵元主要布置在堆芯外圍。燃料柵元的中子價(jià)值隨能量的降低而增大，這是因?yàn)槿剂系牧炎兘孛骐S中子能量的降低而升高。控制棒柵元的中子價(jià)值在能量降低時(shí)迅速下降，原因是控制棒材料對(duì)低能中子的吸收截面很大。水柵元的中子價(jià)值隨能量的降低緩慢降低，原因是水對(duì)低能中子的吸收較強(qiáng)，而快中子能量更高，有更大的可能性經(jīng)慢化后進(jìn)入燃料柵元引起裂變。

圖4 各種柵元在各能群中的中子價(jià)值Fig.4 Neutron importance via energy group in different cells

4 中子價(jià)值對(duì)動(dòng)態(tài)參數(shù)計(jì)算結(jié)果的影響

緩發(fā)中子有效份額βeff和中子代時(shí)間Λ 是反應(yīng)堆動(dòng)態(tài)特性分析中常用到的重要?jiǎng)討B(tài)參數(shù)，它們的定義式中均考慮了中子價(jià)值的影響，βeff和Λ 對(duì)應(yīng)的不考慮中子價(jià)值影響的參數(shù)分別為緩發(fā)中子份額βo 和瞬發(fā)中子壽命lP，本文分析βeff與βo、Λ 與lP存在差異的原因。

4.1 βeff和βo 的比較

緩發(fā)中子有效份額βeff的定義式為：

4.2 Λ 和lP 的比較

中子代時(shí)間Λ 的定義式為：

中子代時(shí)間Λ 的物理意義是：由φ＋（r，E，Ω）作權(quán)重的中子數(shù)和由φ＋（r，E，Ω）作權(quán)重的裂變中子產(chǎn)生率之比。即由φ＋（r，E，Ω）作權(quán)重、每產(chǎn)生一代中子所需的時(shí)間。

瞬發(fā)中子壽命lP的定義式則為：

對(duì)于西安脈沖堆，lP≈93μs，Λ≈36μs。Λ和lP存在明顯的差異，引起這種差異的最重要原因是燃料柵元中的中子價(jià)值較其他柵元中的大。如式（5）所示，Λ 定義式中分子中的中子價(jià)值包括了所有柵元的中子價(jià)值，而分母中的中子價(jià)值則只是燃料柵元中的中子價(jià)值。而lP的定義式中未考慮中子價(jià)值的影響，因此有Λ＜lP。

Λ 和lP的相對(duì)大小，主要是由燃料柵元的平均中子價(jià)值和全堆芯的平均中子價(jià)值的相對(duì)大小決定的，亦即主要是由中子價(jià)值的空間分布決定的。如對(duì)裸堆Godiva和Jezebel［12］，燃料的中子價(jià)值實(shí)際上就是全堆芯的中子價(jià)值，因此二者的Λ 和lP很接近，相對(duì)差別僅為5%左右［13］，而這個(gè)微小的差別是由于中子價(jià)值隨能量的變化引起的，這也同時(shí)說明，只考慮中子價(jià)值隨能量分布的因素，并不會(huì)引起Λ 和lP有很大的差異。

事實(shí)上，由于大多類型的反應(yīng)堆中，燃料柵元的中子價(jià)值均大于或接近全堆芯的平均中子價(jià)值，因此，一般有Λ≤lP。

5 結(jié)論

本文采用柵元程序WIMS 計(jì)算了西安脈沖堆6群群常數(shù)，采用擴(kuò)散程序CITATION 進(jìn)行全堆芯擴(kuò)散計(jì)算，得到了西安脈沖堆的中子價(jià)值在空間和能量上的分布，并獲得以下分析研究結(jié)果：

1）得到了中子價(jià)值隨空間的變化規(guī)律。在選取的15個(gè)柵元中，燃料柵元的中子價(jià)值較高，控制棒和水柵元中的中子價(jià)值相對(duì)較低。燃料柵元在遠(yuǎn)離堆芯或接近控制棒時(shí)，中子價(jià)值會(huì)降低?？熘凶觾r(jià)值在各柵元中的變化較熱中子的價(jià)值變化平緩。

2）得到了各柵元的中子價(jià)值隨能量的變化規(guī)律。燃料柵元的中子價(jià)值隨能量的降低而增大，原因是燃料發(fā)生裂變的可能性隨中子能量的降低而升高?？刂瓢魱旁械闹凶觾r(jià)值隨能量的降低而降低，原因是控制棒材料對(duì)低能中子有很強(qiáng)的吸收作用。水柵元的中子價(jià)值同樣隨能量的降低有緩慢降低，原因是水對(duì)低能中子的吸收較強(qiáng)，高能中子則有更大的可能性經(jīng)慢化后引發(fā)裂變。

3）分析了緩發(fā)中子有效份額βeff和緩發(fā)中子份額βo 存在差異的原因。βeff和βo 的相對(duì)大小，取決于中子價(jià)值隨能量的分布。對(duì)于西安脈沖堆，緩發(fā)中子較瞬發(fā)中子能量低，其平均價(jià)值相對(duì)大，因此有βeff＞βo。

4）分析了中子代時(shí)間Λ 和瞬發(fā)中子壽命lP存在差異的原因。Λ 和lP的相對(duì)大小，主要是由中子價(jià)值隨空間的分布決定的，對(duì)于包含西安脈沖堆在內(nèi)的大多反應(yīng)堆，燃料柵元中的中子價(jià)值較其他柵元中的高，因此有Λ≤lP。

［1］ 杜書華.輸運(yùn)問題的數(shù)值模擬［M］.湖南：湖南科技出版社，1988.

［2］ 謝仲生，鄧力.中子輸運(yùn)理論數(shù)值計(jì)算方法［M］.西安：西北工業(yè)大學(xué)出版社，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.