雙區(qū)次臨界系統(tǒng)的單群Feynman-α 方程的解析解

王子冠，沈 峰

（國核（北京）科學技術研究院有限公司，北京 102209）

Feynman-α方法是中子噪聲分析方法的一種，由Feynman等［1］首先提出，并用來測量增殖系統(tǒng)的反應性。近年來這種方法又被重新研究，并用于ADS 等核能系統(tǒng)次臨界度的在線監(jiān)測［2-3］，也有研究將這種技術應用到核保障中，用于快速檢測貨物中是否存在鈾和钚等其他可裂變元素［4-5］。傳統(tǒng)的Feynman-α方法基于單群點堆模型，且未考慮緩發(fā)中子的影響。但次臨界系統(tǒng)和ADS等核能系統(tǒng)通常由堆芯和反射層構成，并將中子探測器放置于反射層內。對于這類帶有反射層的次臨界系統(tǒng)，基于單群點堆模型的Feynman-α 方法對其物理特性的描述不夠精確?；诖?，本文構建基于雙區(qū)模型的單群Feynman-α方程，并采用構建生成函數的方法推導其解析解，求得中子探測器計數的方差-平均值比值隨探測器計數時間變化的關系。

1 Feynman-α 方法

1.1 基本原理

Feynman-α 方法又被稱為方差-平均值比方法，它通過分析增殖介質內中子探測器計數的方差-平均值比值隨不同探測時間窗長度變化的規(guī)律，計算次臨界系統(tǒng)的α 本征值和反應性等參數。對于一常見的遵循泊松分布的中子源（如Am-Be中子源或252Cf中子源），若將其放置于一不含增殖介質的系統(tǒng)中，則系統(tǒng)內中子探測器計數分布將遵循泊松分布，即中子探測器計數的方差與平均值相同，其方差-平均值比為1。若中子源放置在含有鈾和钚等增殖介質的系統(tǒng)中，則中子將引發(fā)鏈式裂變反應，由1個中子引發(fā)的裂變可能產生多于1個的中子，這將導致系統(tǒng)內的中子之間產生相關性，中子的產生和被探測過程不再是獨立事件，其概率分布也將偏離泊松分布，從而使中子探測器計數的方差-平均值比大于1。這個大于1的部分可用Feynman-Y 方程Y（t）表示：

式中：Z（t）為時間窗長度為t時的中子探測器計數；〈Z（t）〉為探測器計數的一階矩，即探測器計數在探測時間窗長度為t時的平均值；σ（t）2為探測器計數的方差。Y（t）也被稱作Feynman-Y 方程，它描述了中子計數的方差-平均值比值相對于1的偏移量，反映了隨機變量的分布相對于泊松分布的偏移程度。

1.2 中子探測器計數方差和平均值的求解方法

中子產生和被探測屬于隨機過程，因此可運用統(tǒng)計物理學的方法對系統(tǒng)內粒子數量隨時間變化的關系進行研究，求解中子探測器計數的方差和均值，從而得到Feynman-Y 方程。對于隨機過程，利用構造生成函數的方法可求解隨機變量的方差和均值［6］。生成函數是一類以概率密度函數為系數的冪級數，若隨機變量N（t）的概率密度函數為P（N，t），則可定義其生成函數為：

對生成函數求各階導數，可求得隨機變量的各階階乘矩，進而求得隨機變量的期望和方差。例如對于式（2），若先對X 求導，再令X＝1，即可求得隨機變量N（t）一階階乘矩，即N（t）的期望（用尖括號表示）：

若在式（2）中對X 求二階導數，再令X＝1，則可求得隨機變量N（t）二階階乘矩：

則隨機變量的P（N，t）的方差可通過下式算出：

2 雙區(qū)單群Feynman-α方程解析解推導

圖1 雙區(qū)次臨界系統(tǒng)模型示意圖Fig.1 Schematic diagram of two-region subcritical system

構造1個由堆芯和無限厚反射層構成的雙區(qū)次臨界系統(tǒng)模型（圖1），1個遵循泊松分布的中子源位于堆芯內部，源強為S，中子探測器位于反射層內。用P（NA，NB，C，Z，t）表示系統(tǒng)在t時刻處于某種狀態(tài)的概率，即在這種狀態(tài)下，系統(tǒng)中堆芯區(qū)域中子數為NA、反射層區(qū)域中子數為NB、堆內緩發(fā)中子先驅核數為C 且在0到t時間段內探測器中子計數為Z。

假設在單位時間內，反應堆內1個粒子（可能為中子或緩發(fā)中子先驅核）發(fā)生某種反應X的概率（即反應強度）為λX。若系統(tǒng)中已經存在的粒子數為n，則在dt時間內，粒子發(fā)生某種反應X 的概率為：

堆芯區(qū)域（區(qū)域A）中子可能發(fā)生的反應包括：中子被堆內燃料和材料吸收、中子遷移到反射層（區(qū)域B）、中子被核燃料吸收并引發(fā)裂變反應。反射層區(qū)域（區(qū)域B）中子可能發(fā)生的反應有：中子被反射層內材料吸收、中子遷移到堆芯區(qū)域（區(qū)域A）、中子被探測器探測。單位時間內1個粒子發(fā)生各類反應的反應強度用不同字母表示（表1）。

表1 系統(tǒng)內粒子可能發(fā)生的反應類型和反應強度Table 1 Possible reaction types and reaction intensities for different particles in system

在dt時間內，除表1所列的反應類型外，系統(tǒng)內的中子和緩發(fā)中子先驅核還可能均不發(fā)生任何反應，發(fā)生這種情況的概率可寫為：

則在t到dt時間內，系統(tǒng)狀態(tài)的變化可表示為：

將式（8）寫成微分形式，則可構成一描述系統(tǒng)變化情況的前向主方程：

其中：f（k，l）為裂變產生中子數量（k）和緩發(fā)中子先驅核數量（l）的概率密度函數；ps（n）為中子源每次釋放中子數量（n）的概率密度函數。

2.1 求解中子探測器計數期望值表達式

對概率密度函數P（NA，NB，C，Z，t）定義生成函數，有：

由于函數P（NA，NB，C，Z，t）的微分形式相對于其本身更易求得，因此將式（10）寫成微分形式，有：

將式（9）代入式（11）并經過化簡，可得：

為簡化書寫，將式（11）中G（X，Y，V，W，t）簡寫為G。在式（12）中，q（X，V）為概率密度函數f（k，l）的生成函數，定義為：

r（X）為概率密度函數ps（n）的生成函數，定義為：

通過生成函數，可計算出探測器計數的一階階乘矩和二階階乘矩。分別對式（12）中的X、Y、V、W 求導，并令參數X＝Y＝V＝W ＝1，可得系統(tǒng)內粒子數量的一階階乘矩的導數。一階階乘矩的導數描述了系統(tǒng)內各種粒子數的平均值隨時間變化的情況，其中式（15）中的3個微分方程分別描述了系統(tǒng)內堆芯區(qū)域中子數NA、反射層區(qū)域中子數NB、以及緩發(fā)中子先驅核數C 的平均值隨時間變化的情況。式（16）則描述了中子探測器計數Z 的期望值隨時間變化的情況：

其中：

在一源強度不變的次臨界系統(tǒng)中，當時間足夠長時，系統(tǒng)內中子的消耗速率與中子源強度相同，系統(tǒng)達到穩(wěn)態(tài)，系統(tǒng)內各類粒子數的平均值不再隨時間變化，因此各類粒子數的平均值對時間的導數為零，即式（15）中3個微分方程左邊均為0，則式（15）化為了1 個以系統(tǒng)穩(wěn)態(tài)條件下堆芯部分中子數期望值NA、反射層內中子數期望值NB、堆芯內緩發(fā)中子先驅核數期望值C 為變量的三元一次方程組，解此方程組可得：

而通過式（16）可看出，探測器中子計數Z與時間t呈線性關系。對式（16）的時間變量t積分可得中子探測器計數期望值〈Z（t）〉隨探測時間長度變化的表達式：

2.2 求解中子探測器計數方差的表達式

通過生成函數無法直接求解中子探測器計數的方差。但可通過間接方法定義修正二階矩μ 來計算［5］。將修正二階矩μ 定義為：

其中：

當探測時間足夠長，系統(tǒng)處于穩(wěn)態(tài)，與探測器計數無關的各修正二階矩為定值，不隨時間變化，其導數也為零。因此式（23）中6個方程的左邊均為0，則方程組化為以各修正二階矩為變量的六元一次方程組，解此方程可解得與探測器計數無關的各修正二階矩μXX、μXY、μYY、μXV、μYV、μVV。此 六 元 一 次 方 程 組 可 通 過Mathematica等數學軟件求解，但最終得到的解篇幅較長，因此這里未寫出解的具體表達式。

而與探測器計數有關的修正二階矩μXW、μYW、μVW、μWW 對時間的導數也可通過式（22）給出的計算方法得到：

當系統(tǒng)處于穩(wěn)態(tài)時，探測器計數仍隨時間變化而增加，與探測器有關的修正二階矩不為定值，因此式（24）不能轉化為普通方程組求解。為求解式（24），可將其中的3個方程看作是以μXW、μYW、μVW 為變量的微分方程組，利用拉普拉斯變換將之轉化為線性方程組求解。通過這種方法可解得μXW（t）、μYW（t）、μVW（t）的拉普拉斯變換。其中μYW（t）的拉普拉斯變換為：

其中：

H（s）是關于s的三次方程，求解可得到其3個實根，設這3個根分別為-ω1、-ω2、-ω3，即：

對式（26）進行拉普拉斯逆變換，可得μYW（t）：

將式（30）代入式（25），并對t在0到t區(qū)間積分，得：

根據式（1）和式（22）中第1個等式，F(xiàn)eynman-Y 方程可用探測器中子計數的平均值〈Z（t）〉和修正二階矩μWW（t）表示：

因此，將式（21）和式（31）做比值并化簡，即得到雙區(qū)單群Feynman-α方程的解析解：

其中：

3 總結

利用統(tǒng)計物理學方法，通過構造前向主方程，描述了一雙區(qū)次臨界系統(tǒng)的狀態(tài)變化情況，并利用生成函數的性質，求解了系統(tǒng)內中子探測器計數的方差σ2（t）和平均值〈Z（t）〉，進而求得雙區(qū)單群Feynman-α方程的解析解，得到了中子探測器計數的方差-平均值比值隨探測器計數時間變化的關系。目前已有的工作僅限于對雙區(qū)單群Feynman-α 方程的解析解的推導。在下一步的工作中，將嘗試利用蒙特卡羅程序建立計算模型，驗證雙區(qū)Feynman-α方程用于次臨界度求解的有效性，并與單區(qū)單群Feynman-α方法的計算結果進行對比。

［1］ FEYNMAN R P，de HOFFMANN F，SERBER R.Dispersion of the neutron emission in U-235 fission［J］.Journal of Nuclear Energy，1956，3（1）：64-IN10.

［2］ PáZSIT I，KITAMURA Y，WRIGHT J，et al.Calculation of the pulsed Feynman-alpha formulae and their experimental verification［J］.Annals of Nuclear Energy，2005，32（9）：986-1 007.

［3］ TESINSKY M，BERGL?F C，B?CK T，et al.Comparison of calculated and measured reaction rates obtained through foil activation in the subcritical dual spectrum facility YALINA-Booster［J］.Annals of Nuclear Energy，2011，38（6）：1 412-1 417.

［4］ CHERNIKOVA D，PáZSIT I，ZIGUAN W.Application of the two-group-one-region and tworegion-one-group Feynman-alpha formulas in safeguards and accelerator-driven system （ADS）［C］∥Proceeding of ESARDA Meeting 2013.［S.l.］：［s.n.］，2013.

［5］ CHERNIKOVA D，ZIGUAN W，PáZSIT I，et al.A general analytical solution for the varianceto-mean Feynman-alpha formulas for a two-group two-point，a two-group one-point and a onegroup two-point cases［J］.The European Physical Journal Plus，2014，129（11）：1-27.

［6］ PáZSIT I，PáL L.Neutron fluctuations：A treatise on the physics of branching processes［M］.The Netherlands：Elsevier，2007：231-239.