雙區次臨界系統的單群Feynman-α 方程的解析解

王子冠，沈 峰

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

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

1 Feynman-α 方法

1.1 基本原理

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

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

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

中子產生和被探測屬于隨機過程，因此可運用統計物理學的方法對系統內粒子數量隨時間變化的關系進行研究，求解中子探測器計數的方差和均值，從而得到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 雙區單群Feynman-α方程解析解推導

圖1 雙區次臨界系統模型示意圖Fig.1 Schematic diagram of two-region subcritical system

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

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

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

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

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

則在t到dt時間內，系統狀態的變化可表示為：

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

其中：f（k，l）為裂變產生中子數量（k）和緩發中子先驅核數量（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，可得系統內粒子數量的一階階乘矩的導數。一階階乘矩的導數描述了系統內各種粒子數的平均值隨時間變化的情況，其中式（15）中的3個微分方程分別描述了系統內堆芯區域中子數NA、反射層區域中子數NB、以及緩發中子先驅核數C 的平均值隨時間變化的情況。式（16）則描述了中子探測器計數Z 的期望值隨時間變化的情況：

其中：

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

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

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

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

其中：

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

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

當系統處于穩態時，探測器計數仍隨時間變化而增加，與探測器有關的修正二階矩不為定值，因此式（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區間積分，得：

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

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

其中：

3 總結

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