點(diǎn)堆隨機(jī)動(dòng)力學(xué)微分方程的推導(dǎo)

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

(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京100088)

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

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

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

顯然，我們常用的傳統(tǒng)反應(yīng)堆動(dòng)力學(xué)方程組并不具備處理上述物理現(xiàn)象的功能。這是因?yàn)椋绻谙到y(tǒng)中僅有少量的中子、或有隨機(jī)因素時(shí)，則在系統(tǒng)內(nèi)、各空間點(diǎn)上單位體積元內(nèi)中子數(shù)的多少也許不是平均值，而是一個(gè)隨機(jī)事件，系統(tǒng)中的中子數(shù)可能會(huì)非常迅速的增殖，也可能會(huì)全部死絕。

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

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

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

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

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

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

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

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

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

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

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

圖1 田灣核電廠無源啟動(dòng)反應(yīng)性測(cè)量曲線

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

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

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

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

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

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

(a)

(b)

i=2,m+1；

(c)

(d)

n=2,m+1；

(e)

(f)

n=2,m+1；

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

2 點(diǎn)堆隨機(jī)動(dòng)力學(xué)微分方程理論推導(dǎo)

本節(jié)我們將推導(dǎo)和建立引入了隨機(jī)項(xiàng)的點(diǎn)堆隨機(jī)動(dòng)力學(xué)微分方程組。

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

第一種可能的貢獻(xiàn)，俘獲：

(1)

Δt內(nèi)俘獲中子的事件的概率為:

(2)

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

(3)

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

P2=SΔt

(4)

第三種可能的貢獻(xiàn)，Δt內(nèi)裂變釋放ν1個(gè)中子:

(5)

Δt內(nèi)裂變釋放ν1個(gè)中子的概率為：

P3=ΣfυΔty1P(ν1)

(6)

第四種可能的貢獻(xiàn)，Δt內(nèi)裂變釋放ν2個(gè)中子:

(7)

Δt內(nèi)裂變釋放ν2個(gè)中子的概率：

P4=ΣfυΔty1P(ν2)

(8)

第五種可能的貢獻(xiàn)，Δt內(nèi)先驅(qū)核y2衰變：

(9)

Δt內(nèi)先驅(qū)核y2衰變的概率為：

P5=λ2y2Δt

(10)

第六種可能的貢獻(xiàn)，Δt內(nèi)先驅(qū)核y3衰變：

(11)

Δt內(nèi)先驅(qū)核y3衰變的概率為：

P6=λ3y3Δt

(12)

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

則有Δy1、Δy2、Δy3的數(shù)學(xué)期望值為

(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，表達(dá)式如下：

(26)

式中，y(t)=[y1(t),y2(t),…ym+1(t)]是矢量；y1是反應(yīng)堆內(nèi)中子數(shù)；

[y2(t),…ym+1(t)]是相應(yīng)的緩發(fā)中子組隊(duì)先驅(qū)核數(shù)目；

矢量ξ(t)=[ξ1(t),ξ2(t),…ξm+1(t)]由多個(gè)標(biāo)準(zhǔn)的高斯白噪聲組成；

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

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

(27)

(28)

(29)

(30)

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

該方程組被稱為隨機(jī)點(diǎn)堆動(dòng)力學(xué)微分方程，其調(diào)節(jié)部分與點(diǎn)堆動(dòng)力學(xué)方程完全相同。不同的是，式中引入了隨機(jī)項(xiàng)∑gij(y)ξi(t)，或稱為伊藤隨機(jī)項(xiàng)。如果矩陣方程GGT=B=0，則方程組可退化為標(biāo)準(zhǔn)的點(diǎn)堆動(dòng)力學(xué)方程。

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

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

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

(c)式中有bii、yi、i=1三處腳標(biāo)錯(cuò)誤，正確應(yīng)分別為b11、y1、i=2。

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

本文介紹這一物理問題和方程組的推導(dǎo)過程，目的也是想引起業(yè)內(nèi)同行的興趣和關(guān)注，今后，我們將陸續(xù)給出這一方程組的求解方法研究和實(shí)驗(yàn)?zāi)M的結(jié)果。

3 結(jié)論

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