劉建軍,朱家彩，王瑞利，傅學東，任 鍵

(北京應用物理與計算數學研究所,北京100088)

核反應堆點堆動力學方程是一個描述與時間相關的中子數密度和緩發(fā)中子先驅核濃度的常微分方程組。點堆動力學常微分方程組存在單能中子假定、無空間效應等物理描述的局限性。并且該方程組對所有的微觀過程(從核裂變放出中子、介質中的中子輸運、次級核裂變)的每一個環(huán)節(jié)都是離散的隨機過程，最終都要求進行統(tǒng)計處理。

這種處理，只是把本來是隨機過程的中子密度時空分布作為連續(xù)函數的輸運理論、或擴散近似的一種近似描寫[1]。這種處理通常描述的是系統(tǒng)內中子數密度、緩發(fā)中子先驅核濃度、功率狀態(tài)的平均值或期望值。所以，點堆模型核反應動力學方程實際上是確定論的常微分方程，其解也是確定論的平均或期望值。

然而，實際的反應堆動力學過程往往是隨機的，自然，中子密度和緩發(fā)中子先驅核濃度隨著時間變化往往是隨機漲落的。在高功率狀態(tài)，其隨機行為是可以忽略的。但是，在近臨界、低功率狀態(tài)，例如，在弱源條件下的反應堆啟動、臨界安全分析、快中子脈沖堆爆發(fā)脈沖實驗時[2,3]，初始階段將會出現(xiàn)中子密度和緩發(fā)中子增殖過程中的隨機性漲落、以及輻射劑量的不確定度問題，這就特別需要我們從理論和技術上進行把握和掌控，以求準確地評估其范圍和風險。

顯然，我們常用的傳統(tǒng)反應堆動力學方程組并不具備處理上述物理現(xiàn)象的功能。這是因為，如果在系統(tǒng)中僅有少量的中子、或有隨機因素時，則在系統(tǒng)內、各空間點上單位體積元內中子數的多少也許不是平均值，而是一個隨機事件，系統(tǒng)中的中子數可能會非常迅速的增殖，也可能會全部死絕。

雖然，MC方法具有處理此類隨機事件的獨特優(yōu)勢，也最常見。但是，MC方法運算速度慢、人力物力成本相對較高是其主要缺點。

俄羅斯學者在這一問題研究中獨具特色，尤其在反應堆無外源啟動的理論研究和應用實踐中成果顯著，Ю.В.沃爾科夫1988年就提出來隨機動力學的概念。

1992年在他的《弱源條件下的反應堆隨機動力學與核安全研究》[4]一文中簡單介紹了可將隨機微積分應用于描述反應堆物理中的隨機現(xiàn)象，并給出了一組點堆隨機動力學微分方程組。他指出，可以用這一方程組來描述反應堆弱源啟動過程中中子的隨機性漲落物理現(xiàn)象、堆內多群中子持續(xù)時間內的隨機軌跡過程、討論中子數的隨機行為和源強對核安全影響的重要性。

雖然他沒有介紹推導過程，也沒有介紹求解計算方法，我們也難以直觀簡單地判斷其方程的正確性，但是作者創(chuàng)新地在點堆動力學方程中引入了隨機項的建模思路引起了我們極大的關注和興趣。

隨機微積分方法是近年來國際多學科領域處理各種隨機現(xiàn)象的一種熱門的數學方法，日本數學家K.Ito(伊藤清)在這一數學領域的成果引起了國際學術界的極大關注與贊賞。即，求解隨機微分方程最有實用意義的近似方法是伊藤清提出的Ito公式，有了Ito 公式之后，就可以計算一些基本常用的隨機微積分方程。

事實上，這一建模及數學處理方法已經普遍在國際金融、特別是股票交易中漲落情況的預計分析、水利、工程可靠性及不確定性分析等具有隨機事件的領域中都看到諸多成功的研究和應用范例[5-7]。但尚未見到有國內學者在反應堆物理中的應用實踐。

我們借鑒俄羅斯學者Ю.В.沃爾科夫的思想，仔細推導了點堆隨機動力學微分方程組、在點堆動力學方程中引入隨機項的理論研究工作。其目的是，探索快速計算和準確評估反應堆弱源啟動過程中的中子增殖隨機性漲落物理現(xiàn)象及不確定度，并希望引起國內學者的興趣與關注。

本文第一節(jié)，介紹了我們擬研究的物理問題，簡介了Ю.В.沃爾科夫的方程組。第二節(jié)介紹了我們從中子隨機事件的基本假定出發(fā)，推導點堆隨機動力學微分方程組的過程和結果，并指出了Ю.В.沃爾科夫方程組的錯誤。

1 物理問題與沃爾科夫方程組

一般情況下，對于強中子源的核反應堆啟動、運行、或失控時，應用我們熟知的點堆核反應動力學常微分方程組來描述中子的增殖過程是可信的，也不存在隨機問題，只需正確地描述并且考慮最初的無限裂變鏈出現(xiàn)的概率即可。

然而，常用的點堆反應動力學方程在處理和描述含有弱中子源的反應堆啟動、臨界安全分析、脈沖堆爆發(fā)脈沖實驗等物理過程時，就容易出現(xiàn)中子增殖過程中因各種隨機事件發(fā)生隨機性漲落。例如田灣核電站無源(實為燃料棒的自發(fā)裂變源)啟動的反應性測量曲線則直觀地展現(xiàn)了弱源啟動階段反應性與中子增殖過程中的隨機漲落現(xiàn)象[8-10]，如圖1所示。

圖1 田灣核電廠無源啟動反應性測量曲線

此時，對于這種增殖過程存在隨機性漲落，如果采用統(tǒng)計平均的處理方法是存在問題的。

再如，1960年美國學者T.F.Wimett等人公開發(fā)表了在Godiva-Ⅱ脈沖堆進行的弱源爆發(fā)脈沖實驗[3]，也揭示了快堆啟動中的隨機行為和復雜的物理現(xiàn)象。其主要特征是每次爆發(fā)脈沖堆時間呈隨機性，總的爆發(fā)脈沖等待時間概率分布具有規(guī)律性，還會出現(xiàn)提前爆發(fā)脈沖現(xiàn)象，(在脈沖棒未達到預定反應性時就已經提前爆發(fā)脈沖)。

中國工程物理研究院物理與化學研究所也重現(xiàn)了這一現(xiàn)象[11]。

所以，在弱源條件下的核反應堆啟動與安全運行控制時，就要特別關注反應性、或中子數和緩發(fā)中子先驅核的隨機漲落及不確定度，因為，在弱源條件、當k=k(t)時，這種隨機漲落現(xiàn)象及不確定度將影響反應堆的啟動時間，也關系到反應堆的安全風險分析的準確性。

傳統(tǒng)的點堆模型動力學方程實際上是一個數量相互作用的模型系統(tǒng)，特別是中子的數量和緩發(fā)中子先驅核的數量，當物理的動力學系統(tǒng)被標示為數量過程后，在技術上可以把確定論的點堆動力學方程發(fā)展為隨機微分方程系統(tǒng)，即，成為實際的隨機行為過程模型，該隨機微分方程也概括了確定論的點堆動力學方程。

如序言所述，俄羅斯學者Ю.В.沃爾科夫提出了隨機動力學的概念，介紹了可用隨機動力學微積分方程組來描述反應堆物理中的隨機現(xiàn)象。沃爾科夫方程組形式如下：

(a)

(b)

i=2,m+1；

(c)

(d)

n=2,m+1；

(e)

(f)

n=2,m+1；

但是，由于沃爾科夫既沒有介紹推導過程，所以，需要我們通過推導和求證去判斷公式是否正確。

2 點堆隨機動力學微分方程理論推導

本節(jié)我們將推導和建立引入了隨機項的點堆隨機動力學微分方程組。

考慮一個十分小的時間間隔Δt，Δt時間內有2+Q種不同事件單獨發(fā)生。因為Δt足夠小，可以假定任何兩種事件不同時發(fā)生，y1(t)表示系統(tǒng)內中子數隨時間的變化，考慮存在兩組緩發(fā)中子β2、β3，則Δt內有6種不同事件單獨發(fā)生，

第一種可能的貢獻，俘獲：

(1)

Δt內俘獲中子的事件的概率為:

(2)

第二種可能的貢獻，定常中子源S在Δt內發(fā)源射中子:

(3)

定常中子源S在Δt內發(fā)射的概率為:

P2=SΔt

(4)

第三種可能的貢獻，Δt內裂變釋放ν1個中子:

(5)

Δt內裂變釋放ν1個中子的概率為：

P3=ΣfυΔty1P(ν1)

(6)

第四種可能的貢獻，Δt內裂變釋放ν2個中子:

(7)

Δt內裂變釋放ν2個中子的概率：

P4=ΣfυΔty1P(ν2)

(8)

第五種可能的貢獻，Δt內先驅核y2衰變：

(9)

Δt內先驅核y2衰變的概率為：

P5=λ2y2Δt

(10)

第六種可能的貢獻，Δt內先驅核y3衰變：

(11)

Δt內先驅核y3衰變的概率為：

P6=λ3y3Δt

(12)

上述式中，y1(t)表示系統(tǒng)內中子數隨時間的變化；先驅核的份額分別為β2、β3；裂變產生次級中子的數目ν1和ν2的分布為P(ν1)和P(ν2)。y2、y3分別為緩發(fā)中子先驅核的數目；L是中子壽命；Pi為產生事件i的概率；υ為中子速度；Σa為中子的宏觀吸收截面；Σf為中子的宏觀裂變截面；Σγ為中子的宏觀俘獲截面。

則有Δy1、Δy2、Δy3的數學期望值為

(13)

則有：

(14)

平均協(xié)方差

(15)

λ2y2Δt+λ3y3Δt+ΣfvΔty1P(ν2)[-1+(1-β)ν2]2

λ2y2Δt+λ3y3Δt

(16)

B12Δt=P1[Δy1Δy2]1+P2[Δy1Δy2]2+…+P6[Δy1Δy2]6

=0+0+ΣfvΔty1P(ν1)[-1+(1-β)ν1]β2ν1

+ΣfvΔty1P(ν2)[-1+(1-β)ν2]β3ν2-λ2y2Δt

(17)

(18)

(19)

比較可知，

B21=B12

(20)

B31=B13

(21)

(22)

(23)

B23Δt=P1[Δy2Δy3]1+P2[Δy2Δy3]2+…+P6[Δy2Δy3]6

(24)

同樣有：

B32=B23

(25)

以此類推，我們就不難得到與沃爾科夫方程相似的方程組及其元素Bij，表達式如下：

(26)

式中，y(t)=[y1(t),y2(t),…ym+1(t)]是矢量；y1是反應堆內中子數；

[y2(t),…ym+1(t)]是相應的緩發(fā)中子組隊先驅核數目；

矢量ξ(t)=[ξ1(t),ξ2(t),…ξm+1(t)]由多個標準的高斯白噪聲組成；

式(26)中右側部分白噪聲的存在表示中子出生和死亡等基本過程的特征，其特征時間?L。

gij是解矩陣方程GGT=B的元素，維數為(m+1)(m+1)的矩陣，Bij的元素為：

(27)

(28)

(29)

(30)

到此，我們得到了完整的隨機中子微分方程組式(26)～方程組式(30)。

該方程組被稱為隨機點堆動力學微分方程，其調節(jié)部分與點堆動力學方程完全相同。不同的是，式中引入了隨機項∑gij(y)ξi(t)，或稱為伊藤隨機項。如果矩陣方程GGT=B=0，則方程組可退化為標準的點堆動力學方程。

方程組可以描述反應堆內有少量的源中子和存在多群緩發(fā)中子先驅核的連續(xù)時間內的隨機過程的軌跡。為了評估反應堆的安全性，或發(fā)生事故時的輻射劑量，則必須弄清y(t)=[y1(t),y2(t),…ym+1(t)]的時間過程，

還必須解決中間概率下的多為擴散過程y(t)的邊界問題，該問題都可通過上述點堆隨機動力學微分方程組來描述。在引進了隨機項后，就可以來描述中子增殖過程中的隨機漲落現(xiàn)象，而不僅僅是期望值。

我們把本文得到的方程組與沃爾科夫的方程組和元素Bij進行了仔細的對比。我們發(fā)現(xiàn)不知是沃爾科夫筆誤還是印刷錯誤，在他的方程組與腳標存在以下多處錯誤：

(c)式中有bii、yi、i=1三處腳標錯誤，正確應分別為b11、y1、i=2。

通過上述工作，我們發(fā)現(xiàn)了沃爾科夫點堆隨機動力學方程組的某些錯誤，得到了正確的方程組，這將為我們后續(xù)對方程組的求解方法研究和實驗的模擬奠定了可靠的基礎。我們研究認為，雖然沃爾科夫所給的方程組存在一些筆誤或錯誤，但他的物理思想是科學創(chuàng)新的。而且，隨機反應堆動力學方程也將是今后反應堆物理中值得研究和發(fā)展的方向之一。

本文介紹這一物理問題和方程組的推導過程，目的也是想引起業(yè)內同行的興趣和關注，今后，我們將陸續(xù)給出這一方程組的求解方法研究和實驗模擬的結果。

3 結論

針對反應堆物理中的隨機問題的研究，通過調研和借鑒國外學者的理念，我們嘗試著將隨機微積分理論引入于反應堆物理，在點堆模型假定下，仔細推導了可描述隨機現(xiàn)象的點堆隨機動力學微分方程組，發(fā)現(xiàn)并糾正了沃爾科夫方程組的錯誤。這將為我們進行方程組的求解方法研究和實驗的模擬奠定了可靠的基礎。