王 喆 洪振英

（北京應用物理與計算數(shù)學研究所 北京 100094）

中子引發(fā)裂變鏈概率的演化過程模擬

王 喆 洪振英

（北京應用物理與計算數(shù)學研究所 北京 100094）

在相對速度空間建立中子引發(fā)裂變鏈概率所滿足的與時間相關的微分-積分方程，基于多群SN方法開發(fā)動態(tài)數(shù)值程序(Dynamic Segment Number Probability, DSNP)，分析了動態(tài)計算的收斂性，并對動態(tài)系統(tǒng)的裂變鏈概率演化過程進行數(shù)值模擬。模擬計算表明，DSNP程序與Partisn程序的計算結(jié)果均一致；臨界狀態(tài)附近存在大量的有限裂變鏈，使得引發(fā)概率的動態(tài)演化結(jié)果高于穩(wěn)態(tài)計算結(jié)果，在Baker動態(tài)流場模型上，第一臨界點后1 μs的范圍內(nèi)計算結(jié)果最大差異約為300%。隨著裂變系統(tǒng)反應性增加，有限裂變鏈的貢獻逐漸減弱，持續(xù)裂變鏈占優(yōu)，引發(fā)概率的動態(tài)演化曲線與穩(wěn)態(tài)結(jié)果逐漸重合，差別小于5%，表明系統(tǒng)中子引發(fā)自持裂變的能力趨于穩(wěn)定，此時動態(tài)引發(fā)概率的時間積分結(jié)果比穩(wěn)態(tài)結(jié)果高5%-35%。高濃鈾模型上的數(shù)值模擬驗證了DSNP程序的準確性，該程序可定量計算動態(tài)系統(tǒng)的引發(fā)概率，相對于穩(wěn)態(tài)方法，DSNP程序能夠更為準確地描述裂變系統(tǒng)點火概率的演化過程。

持續(xù)裂變鏈概率，動態(tài)裂變系統(tǒng)，多群SN方法

中子引發(fā)的持續(xù)裂變鏈概率是含裂變材料增殖系統(tǒng)的固有參量，是臨界安全分析、反應堆的啟動、脈沖堆爆發(fā)脈沖等待時間概率分布等工作中均需計算的物理量。描述該量的控制方程早在20世紀60年代由Bell等[1-2]建立，當時間趨于無窮時，持續(xù)裂變鏈概率滿足與時間無關的微分-積分方程，該穩(wěn)態(tài)方程及其解僅適用于描述固定力學狀態(tài)下的持續(xù)裂變鏈概率，反映裂變鏈在無限長時間內(nèi)的行為。

對于力學狀態(tài)隨時間變化的裂變系統(tǒng)而言，中子與原子核相互作用的隨機性決定了實際產(chǎn)生的任何裂變鏈僅能存在有限時間，中子引發(fā)的裂變鏈概率與增殖系統(tǒng)內(nèi)部中子場相類似，表現(xiàn)為一種隨機演化的過程，需要由與時間有關的微分-積分方程來描述，并利用末態(tài)條件進行求解。國內(nèi)外的研究工作表明[3-7]，一直以來，這個領域的研究工作大多局限于定態(tài)構型的穩(wěn)態(tài)方程求解上，國內(nèi)早在20世紀70年代末開發(fā)了定態(tài)SN程序（計算結(jié)果未公開發(fā)表）。國際上在1976年編制了SNP程序[3]，后來陸續(xù)發(fā)展了多個程序來進行確定論求解。如Partisn[5]、Ardra[6]和Amtran程序[7]，以及基于概率論的蒙卡程序Mercury[8]；直到2009年，含時方程的求解才見報道，Baker[9]在其文章中介紹了動態(tài)裂變系統(tǒng)中的裂變鏈概率演化計算結(jié)果。

裂變鏈初期發(fā)展具有顯著的統(tǒng)計漲落性質(zhì)，只有少數(shù)裂變鏈能夠持續(xù)發(fā)展下去。為了進一步研究動態(tài)增殖系統(tǒng)中裂變鏈概率的隨機演化過程，分析有限裂變鏈對引發(fā)概率的貢獻，本文在相對速度空間建立中子引發(fā)概率所滿足的主方程，基于多群SN方法自主開發(fā)動態(tài)數(shù)值程序(Dynamic Segment Number Probability, DSNP)，在一維球?qū)ΨQ坐標系下對中子引發(fā)裂變鏈概率的主方程進行求解，初步討論了有限裂變鏈對裂變鏈概率動態(tài)演化結(jié)果的影響程度。由于動態(tài)系統(tǒng)裂變鏈概率演化過程的數(shù)值求解尚屬首次，首先在定態(tài)構型上驗證動態(tài)計算收斂的正確性，其次研究動態(tài)系統(tǒng)的裂變鏈概率演化與穩(wěn)態(tài)計算結(jié)果的差異，分析臨界狀態(tài)附近有限裂變鏈對引發(fā)概率的貢獻。數(shù)值驗證所用模型取自文獻[9]中的定態(tài)測試幾何結(jié)構和動態(tài)假想流場，結(jié)果表明，在能量分群、中子離散方向、各項異性截斷與迭代誤差判據(jù)各自不同的條件下，DSNP程序的計算結(jié)果與國外確定論方法[8-9]給出的結(jié)果是一致的，臨界狀態(tài)附近大量有限裂變鏈的存在，使引發(fā)概率的動態(tài)演化結(jié)果高于穩(wěn)態(tài)計算結(jié)果。隨著裂變系統(tǒng)反應性的增加，有限裂變鏈的貢獻逐漸減弱，持續(xù)裂變鏈占優(yōu)，引發(fā)概率的動態(tài)演化曲線與穩(wěn)態(tài)結(jié)果逐漸重合，系統(tǒng)中子引發(fā)自持裂變的能力趨于穩(wěn)定。模擬結(jié)果驗證了DSNP程序的準確性，可以用來定量計算動態(tài)系統(tǒng)的引發(fā)概率，與穩(wěn)態(tài)方法相比，DSNP程序提高了對裂變系統(tǒng)的點火概率問題的模擬和論證能力，為實際工程應用提供了更為準確的設計與評估手段。

1 與時間有關的微分-積分方程

t時刻初始中子未引發(fā)的持續(xù)裂變鏈概率等于t′時刻的后代中子未引發(fā)的持續(xù)裂變鏈概率，這是建立微分-積分方程的基本守恒關系。中子與核的相互作用分別用ΣT、Σf、Σs、Σc表示總截面、裂變截面、散射截面和俘獲截面描述，則守恒關系表示為：

等號右邊對裂變項的二項式作展開，并利用截面關系式ΣT=Σf+Σs+Σc，得到方程：

一維球?qū)ΨQ坐標系下的多群方程可以表示為：

方程滿足的定解條件為

t=tf時，ωg(r, μ, tf)=ε (0≤ε≤1)；

μ≥0時，ωg(R, μ, t)=0；

交界面處，ωg(r+, μ, t)= ωg(r-, μ, t)；

中心軸對稱處，ωg(0, μ, t)= ωg(0, μ', t)，μ與μ'滿足μ=-μ'，且

2 數(shù)值計算方法

離散縱標(SN)法是研究粒子輸運問題及求解相關方程的有效數(shù)值計算方法之一。本文選擇高斯-勒讓德求積組，對空間變量r、角方向μ進行離散，將r-μ平面分割為網(wǎng)格形式，Δk+1/2,m=(rk, rk+1, μm-1/2, μm+1/2)，k=0, 1, 2, …, K-1，m=1, 2, …, N；時間也采用類似形式表示。空間變量r、角方向μ、時間變量t的中心量與邊界量的符號表示見表1。r1=0為系統(tǒng)中心點，R為系統(tǒng)外邊界，tn+1對應于實際力學狀態(tài)的時刻，時間間隔Δtn+1=(tn+3/2-tn+1/2)，n=0, 1, 2, …, N-1。

方程(4)右邊源項用符號Qg表示，在空間、角方向和時間上積分，得到主方程的差分方程形式：

表1 中心量與邊界量符號Table 1 Signs of center and boundary variable.

3 數(shù)值驗證

主方程(4)按照時間反向進行求解，其初始條件實際為末態(tài)條件。對于一個末態(tài)為超臨界的系統(tǒng)，末態(tài)上的穩(wěn)態(tài)解可以作為末態(tài)條件，也可由定義，以tf時刻的一個中子在tf時刻的存活概率“1”作為末態(tài)條件。

3.1 動態(tài)計算的收斂性驗證

系統(tǒng)的力學狀態(tài)不發(fā)生變化時，無論初值如何選取，只要計算時間足夠長（視問題而定），求解含時方程得到的裂變鏈概率演化過程，其漸進解與穩(wěn)態(tài)解相同。驗證計算在高濃鈾裸球模型上進行，模型的質(zhì)量密度與尺寸及本征值數(shù)據(jù)[9]見表2，同位素原子比U235為9.38×10-1，U238為6.20×10-2。

利用DSNP程序計算中子持續(xù)裂變鏈概率的時間、空間、能量、中子角方向的聯(lián)合分布，再以歸一化的源分布為權重計算裂變鏈概率的加權平均值，即一個中子引發(fā)的點火概率PS(t)，計算公式如下：

源分布中裂變譜為Watt譜，空間均勻分布(rsource=rsphere)，歸一化條件為

動態(tài)計算的時間為0-40 sh（sh為時間單位，1sh=0.01 μs），t=40 sh時，末態(tài)條件分別為穩(wěn)態(tài)方程的解（settle形式）和常分布“1”（given形式）。

表2 高濃鈾裸球模型Table 2 Highly enriched uranium bare spheres.

文獻[9]給出了LANL的SN輸運程序Partisn的計算結(jié)果，能量分群為30群，中子離散角方向數(shù)為S16，勒讓德展開取P4，裂變釋放中子數(shù)最大值取8；DSNP程序計算時，能量分群為24群，中子離散角方向數(shù)為S16，勒讓德展開取P3，裂變釋放中子數(shù)最大值取10。計算結(jié)果如圖1所示。

圖1 Partisn程序和DSNP程序計算的高濃鈾裸球模型(a) 中子引發(fā)概率的穩(wěn)態(tài)結(jié)果，(b) 動態(tài)結(jié)果Fig.1 Static results (a) and dynamic results (b) of probability of initiation on highly-enriched-uranium bare sphere by Partisn code and DSNP code.

圖1 (a)中六條水平方向的線表示初值為穩(wěn)態(tài)方程的解，由于末態(tài)條件本身就是穩(wěn)態(tài)方程的解，在構型不變的情況下，整個計算過程中裂變鏈概率不再隨時間變化；末態(tài)條件為“1”時，隨著計算向t=0方向進行，非平衡的末態(tài)條件導致的瞬變過程逐漸消失，動態(tài)引發(fā)概率逐漸趨于其穩(wěn)態(tài)解，如表3所示。

表3 t=0 sh時的一個中子引發(fā)概率Table 3 Single particle initiation probabilities at t=0 sh.

計算結(jié)果表明，無論初值如何給定，定態(tài)構型上中子裂變鏈概率的動態(tài)演化漸進解與穩(wěn)態(tài)解是相等的；DSNP程序與Partisn程序的計算結(jié)果非常接近，差別小于2%，在系統(tǒng)反應性較大時，計算結(jié)果的差異不到1%。面密度為165 g?cm-2的構型是次臨界的，時間趨于無窮時中子存活概率為0。

3.2 動態(tài)流場的概率演化計算

據(jù)文獻[9]，動態(tài)假想流場為半徑17.25 cm的球模型，材料為U235/U238，密度為15.0 g?cm-3。同位素核所占質(zhì)量百分比隨時間變化的關系如下：

式中，CMOD是一個隨時間變化的分段函數(shù)，使系統(tǒng)反應性隨時間變化，以此來表征動態(tài)系統(tǒng)，如圖2所示。

圖2 動態(tài)系統(tǒng)的組份與臨界度參數(shù)Fig.2 Concentration/criticality variation on dynamic system.

在這個反應性隨時間變化的動態(tài)系統(tǒng)中，計算每個固定力學狀態(tài)上的中子持續(xù)裂變鏈概率以及裂變鏈概率的動態(tài)演化分布，再以歸一化的源分布為權重計算裂變鏈概率的加權平均值，即一個中子引發(fā)的點火概率PS(t)。動態(tài)計算時末態(tài)條件分別為穩(wěn)態(tài)方程的解（settle形式）和常分布“1”（given形式）。以穩(wěn)態(tài)解為末態(tài)條件時的起始計算時刻為第二臨界點附近，以常分布“1”為末態(tài)條件時的起始計算時刻分別為800 sh、1000 sh。

Partisn程序的能量分群為21群、中子離散角方向數(shù)為S20、勒讓德展開取P3、裂變釋放中子數(shù)最大值取6，收斂判據(jù)為10-6[9]；DSNP程序的能量分群為24群、中子離散角方向數(shù)為S16、勒讓德展開取P3、裂變釋放中子數(shù)最大值取10，收斂判據(jù)為10-8；計算結(jié)果如圖3所示，兩個程序的計算結(jié)果具有較好的一致性。

圖3 Partisn程序和DSNP程序計算動態(tài)系統(tǒng)的中子引發(fā)概率結(jié)果(a)和t=800 sh的動態(tài)結(jié)果比較圖(b)Fig.3 Probability of initiation on varying concentration (a) and comparison of dynamic results at t=800 sh (b) by Partisn code and DSNP code.

3.3 結(jié)果分析與討論

在文獻[9]中使用的動態(tài)假想流場是迄今唯一公開報道的“基準”動態(tài)模型，本文對該算例的計算結(jié)果與其高度一致，對裂變鏈概率演化過程的物理特征分析與文獻[9]也基本一致（圖4）。裂變鏈概率的動態(tài)演化過程表明臨界狀態(tài)附近有限裂變鏈使得引發(fā)概率的動態(tài)演化結(jié)果高于穩(wěn)態(tài)計算結(jié)果。

圖4 動態(tài)與穩(wěn)態(tài)計算結(jié)果的差異比較(a) 中子引發(fā)概率，(b) 卷積結(jié)果Fig.4 Difference of static and dynamic results. (a) Probability of initiation, (b) Cumulative probability

本算例所采用的動態(tài)假想流場是關于t=450 sh對稱的，穩(wěn)態(tài)計算給出的點火概率曲線也是對稱的，因為穩(wěn)態(tài)方程描述的是時間趨于無窮時的中子持續(xù)裂變鏈概率，每個力學狀態(tài)下的計算結(jié)果同時間序列無關；第一臨界點前穩(wěn)態(tài)計算結(jié)果均為“0”，表明處于次臨界狀態(tài)的系統(tǒng)不會產(chǎn)生持續(xù)裂變鏈，而動態(tài)演化的結(jié)果不為“0”，表明次臨界狀態(tài)下的系統(tǒng)中存在大量的有限裂變鏈，可能存活一定的時間，少量甚至存活至第二臨界點以后[9]，第一臨界點后至少約1 μs的時間范圍內(nèi)，動態(tài)結(jié)果也比穩(wěn)態(tài)結(jié)果高，見圖4(a)，二者的曲線點是分離的，最大差異約300%；總的來說，臨界狀態(tài)附近（無論是臨界前還是臨界后）動態(tài)演化結(jié)果均高于穩(wěn)態(tài)計算結(jié)果，原因在于臨界狀態(tài)附近存在大量的有限裂變鏈，既形成了中子數(shù)目的統(tǒng)計漲落現(xiàn)象，又因?qū)χ凶右l(fā)概率的貢獻導致了動態(tài)演化結(jié)果高于穩(wěn)態(tài)計算結(jié)果。

由于裂變鏈初期具有顯著的隨機漲落性質(zhì)，只有部分裂變鏈能夠形成持續(xù)裂變鏈發(fā)展下去，其它裂變鏈在達到最大鏈長以前將會中止，中子引發(fā)概率主要來自于持續(xù)裂變鏈，也就是說，如果系統(tǒng)已經(jīng)度過了裂變鏈增殖的早期階段，系統(tǒng)中的有限裂變鏈數(shù)目減少，持續(xù)裂變鏈逐漸占優(yōu)，中子引發(fā)概率的動態(tài)演化結(jié)果應該等于或接近穩(wěn)態(tài)方程給出的持續(xù)裂變鏈概率結(jié)果；如圖4(a)所示，高超臨界狀態(tài)下，由于裂變鏈增殖時間仍然足夠長，中子引發(fā)概率的動態(tài)結(jié)果逼近于穩(wěn)態(tài)結(jié)果，動態(tài)演化曲線與穩(wěn)態(tài)結(jié)果逐漸重合，二者最大差異低于5%。

由于次臨界狀態(tài)下的穩(wěn)態(tài)計算結(jié)果為0，因此即使在臨界附近動態(tài)卷積結(jié)果較?。s10-2量級），仍遠高于穩(wěn)態(tài)的卷積結(jié)果[9]，隨著系統(tǒng)反應性的增加，動態(tài)卷積結(jié)果與穩(wěn)態(tài)結(jié)果的差異逐漸減少，如圖4(b)所示，在300 sh以后，二者差異為5%-35%。

早期引發(fā)的裂變鏈概率的時間積分效應，其動態(tài)演化結(jié)果與穩(wěn)態(tài)結(jié)果的差異在實際工程應用中是不能忽略的，該量的大小對于評估裂變系統(tǒng)的點火概率、安全性以及確定外中子源技術要求都有非常敏感而且重要的作用。

4 結(jié)語

本文在相對速度空間建立了主方程，基于多群SN方法開發(fā)了DSNP程序，首次實現(xiàn)了對中子引發(fā)裂變鏈概率含時主方程的求解。通過對動態(tài)計算收斂性和裂變鏈概率演化過程所進行的數(shù)值模擬，得到以下結(jié)論：

(1) 自主開發(fā)的DSNP程序與國外確定論方法的計算結(jié)果符合；

(2) 裂變鏈概率動態(tài)演化過程表明，由于有限裂變鏈的存在，次臨界狀態(tài)下的動態(tài)演化結(jié)果不為0，早期的中子引發(fā)概率動態(tài)演化結(jié)果顯著高于穩(wěn)態(tài)結(jié)果，隨著系統(tǒng)反應性增加，動態(tài)結(jié)果逐漸接近于穩(wěn)態(tài)計算結(jié)果；

(3) 中子引發(fā)概率的動態(tài)演化結(jié)果在臨界狀態(tài)附近約為穩(wěn)態(tài)計算結(jié)果的3倍，高超臨界時二者的差異約為5%；引發(fā)概率的時間卷積結(jié)果在臨界附近的差別很大，可能達到幾十倍甚至一百倍，隨著系統(tǒng)反應性的增加，差異逐漸減小，高超臨界態(tài)下動態(tài)與穩(wěn)態(tài)結(jié)果的差異為5%-35%；早期引發(fā)概率的時間積分結(jié)果、動態(tài)演化結(jié)果與穩(wěn)態(tài)結(jié)果的差異不可忽略。

1 Hansen G E. Assembly of fissionable material in the presence of weak neutron source[J]. Nuclear Science and Engineering, 1960, 8: 709-719

2 Bell G I. On the stochastic theory of neutron transport[J]. Nuclear Science and Engineering, 1965, 21: 390-401

3 Bell G I, Lee C E. On the probability of initiating a persistent fission chain[R]. Los Alamos National Laboratory, 1976

4 劉建軍, 王喆, 張本愛. 隨機中子輸運方程數(shù)值解與近似解的對比分析[J]. 核科學與工程, 2006, 26(2): 118-121 LIU Jianjun, WANG Zhe, ZHANG Ben’ai. Research of stochastic neutron transport equation’s numerial and analytic solution[J]. Chinese Journal of Nuclear Science and Engineering, 2006, 26(2):118-121

5 Alcouffe R E, Baker R S, Turner S A, et al. PARTISN: a time-dependent, parallel neutral particle transport code system[R]. Los Alamos National Laboratory, 2002

6 Hannebutte U, Brown P N. Ardra: scalable parallel code system to perform neutron and radiation transport calculations[R]. Lawrence Livermore National Laboratory, 1999

7 Clouse C J. Parallel deterministic neutron transport with AMR[A]. In: Graziani F R, Springer-Verlag. Computational Methods in Transport[M]. 2006: 499-512

8 Greenman G M, Procassini R J, Clouse C J. A Monte Carlo method for calculating initiation probability[C]. Joint International Topical Meeting on Mathematics, Computation, and Supercomputing in Nuclear Applications, Monterrey, CA, April 15-19, 2007

9 Baker R S. Deterministic methods for time-dependent stochastic neutron transport[C]. International Conference on Mathematics, Computational Methods & Reactor Physics (M&C 2009), Saratoga Springs, New York, May 3-7, 2009

CLCO571.43+7, TL329

Evolvement simulation of the probability of neutron-initiating persistent fission chain

WANG Zhe HONG Zhenying
(Beijing Institute of Applied Physics and Computational Mathematics, Beijing 100094, China)

Background:Probability of neutron-initiating persistent fission chain, which has to be calculated in analysis of critical safety, start-up of reactor, burst waiting time on pulse reactor, bursting time on pulse reactor, etc., is an inherent parameter in a multiplying assembly.Purpose:We aim to derive time-dependent integro-differential equation for such probability in relative velocity space according to the probability conservation, and develop the deterministic code Dynamic Segment Number Probability (DSNP) based on the multi-group SNmethod.Methods:The reliable convergence of dynamic calculation was analyzed and numerical simulation of the evolvement process of dynamic probability for varying concentration was performed under different initial conditions.Results:On Highly Enriched Uranium (HEU) Bare Spheres, when the time is long enough, the results of dynamic calculation approach to those of static calculation. The most difference of such results between DSNP and Partisn code is less than 2%. On Baker model, over the range of about 1 μs after the first criticality, the most difference between the dynamic and static calculation is about 300%. As for a super critical system, the finite fission chains decrease and the persistent fission chains increase as the reactivity aggrandizes, the dynamic evolvement curve of initiation probability is close to the static curve within the difference of 5% when the Keffis more than 1.2. The cumulative probability curve also indicates that the difference of integral results between the dynamic calculation and the static calculation decreases from 35% to 5% as the Keffincreases. This demonstrated that the ability of initiating a self-sustaining fission chain reaction approaches stabilization, while the former difference (35%) showed the important difference of the dynamic results near the first criticality with the static ones. The DSNP code agrees well with Partisn code.Conclusions:There are large numbers of finite fission chains near the first criticality, which can survive until after the second criticality. So the results of dynamic calculation here will be greater than those of static calculation. The numerical simulation on HEU Bare Spheres and Baker model validated the accuracy of the DSNP code, which can calculate initiation probability of dynamic system. Relative to the static method, the DSNP code can describe perfectly the evolvement process of probability of ignition of fissile system.

Probability of persistent fission chain, Dynamic fissile system, Multi-group SNmethod

O571.43+7，TL329

10.11889/j.0253-3219.2014.hjs.37.050602

中國工程物理研究院科學技術發(fā)展基金項目(No.2013B0103017)資助

王喆，男，1973年4月出生，1999年于吉林大學獲碩士學位，副研究員，主要從事核反應堆物理設計工作

2014-01-06，

2014-01-28