王子冠，沈 峰

（國(guó)核（北京）科學(xué)技術(shù)研究院有限公司，北京 102209）

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

1 Feynman-α 方法

1.1 基本原理

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

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

1.2 中子探測(cè)器計(jì)數(shù)方差和平均值的求解方法

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

對(duì)生成函數(shù)求各階導(dǎo)數(shù)，可求得隨機(jī)變量的各階階乘矩，進(jìn)而求得隨機(jī)變量的期望和方差。例如對(duì)于式（2），若先對(duì)X 求導(dǎo)，再令X＝1，即可求得隨機(jī)變量N（t）一階階乘矩，即N（t）的期望（用尖括號(hào)表示）：

若在式（2）中對(duì)X 求二階導(dǎo)數(shù)，再令X＝1，則可求得隨機(jī)變量N（t）二階階乘矩：

則隨機(jī)變量的P（N，t）的方差可通過(guò)下式算出：

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

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

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

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

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

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

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

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

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

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

2.1 求解中子探測(cè)器計(jì)數(shù)期望值表達(dá)式

對(duì)概率密度函數(shù)P（NA，NB，C，Z，t）定義生成函數(shù)，有：

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

將式（9）代入式（11）并經(jīng)過(guò)化簡(jiǎn)，可得：

為簡(jiǎn)化書(shū)寫(xiě)，將式（11）中G（X，Y，V，W，t）簡(jiǎn)寫(xiě)為G。在式（12）中，q（X，V）為概率密度函數(shù)f（k，l）的生成函數(shù)，定義為：

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

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

其中：

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

而通過(guò)式（16）可看出，探測(cè)器中子計(jì)數(shù)Z與時(shí)間t呈線性關(guān)系。對(duì)式（16）的時(shí)間變量t積分可得中子探測(cè)器計(jì)數(shù)期望值〈Z（t）〉隨探測(cè)時(shí)間長(zhǎng)度變化的表達(dá)式：

2.2 求解中子探測(cè)器計(jì)數(shù)方差的表達(dá)式

通過(guò)生成函數(shù)無(wú)法直接求解中子探測(cè)器計(jì)數(shù)的方差。但可通過(guò)間接方法定義修正二階矩μ 來(lái)計(jì)算［5］。將修正二階矩μ 定義為：

其中：

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

而與探測(cè)器計(jì)數(shù)有關(guān)的修正二階矩μXW、μYW、μVW、μWW 對(duì)時(shí)間的導(dǎo)數(shù)也可通過(guò)式（22）給出的計(jì)算方法得到：

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

其中：

H（s）是關(guān)于s的三次方程，求解可得到其3個(gè)實(shí)根，設(shè)這3個(gè)根分別為-ω1、-ω2、-ω3，即：

對(duì)式（26）進(jìn)行拉普拉斯逆變換，可得μYW（t）：

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

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

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

其中：

3 總結(jié)

利用統(tǒng)計(jì)物理學(xué)方法，通過(guò)構(gòu)造前向主方程，描述了一雙區(qū)次臨界系統(tǒng)的狀態(tài)變化情況，并利用生成函數(shù)的性質(zhì)，求解了系統(tǒng)內(nèi)中子探測(cè)器計(jì)數(shù)的方差σ2（t）和平均值〈Z（t）〉，進(jìn)而求得雙區(qū)單群Feynman-α方程的解析解，得到了中子探測(cè)器計(jì)數(shù)的方差-平均值比值隨探測(cè)器計(jì)數(shù)時(shí)間變化的關(guān)系。目前已有的工作僅限于對(duì)雙區(qū)單群Feynman-α 方程的解析解的推導(dǎo)。在下一步的工作中，將嘗試?yán)妹商乜_程序建立計(jì)算模型，驗(yàn)證雙區(qū)Feynman-α方程用于次臨界度求解的有效性，并與單區(qū)單群Feynman-α方法的計(jì)算結(jié)果進(jìn)行對(duì)比。

［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.