劉建軍,朱家彩,王瑞利,傅學(xué)東,任 鍵
(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京100088)
核反應(yīng)堆點堆動力學(xué)方程是一個描述與時間相關(guān)的中子數(shù)密度和緩發(fā)中子先驅(qū)核濃度的常微分方程組。點堆動力學(xué)常微分方程組存在單能中子假定、無空間效應(yīng)等物理描述的局限性。并且該方程組對所有的微觀過程(從核裂變放出中子、介質(zhì)中的中子輸運、次級核裂變)的每一個環(huán)節(jié)都是離散的隨機過程,最終都要求進行統(tǒng)計處理。
這種處理,只是把本來是隨機過程的中子密度時空分布作為連續(xù)函數(shù)的輸運理論、或擴散近似的一種近似描寫[1]。這種處理通常描述的是系統(tǒng)內(nèi)中子數(shù)密度、緩發(fā)中子先驅(qū)核濃度、功率狀態(tài)的平均值或期望值。所以,點堆模型核反應(yīng)動力學(xué)方程實際上是確定論的常微分方程,其解也是確定論的平均或期望值。
然而,實際的反應(yīng)堆動力學(xué)過程往往是隨機的,自然,中子密度和緩發(fā)中子先驅(qū)核濃度隨著時間變化往往是隨機漲落的。在高功率狀態(tài),其隨機行為是可以忽略的。但是,在近臨界、低功率狀態(tài),例如,在弱源條件下的反應(yīng)堆啟動、臨界安全分析、快中子脈沖堆爆發(fā)脈沖實驗時[2,3],初始階段將會出現(xiàn)中子密度和緩發(fā)中子增殖過程中的隨機性漲落、以及輻射劑量的不確定度問題,這就特別需要我們從理論和技術(shù)上進行把握和掌控,以求準確地評估其范圍和風(fēng)險。
顯然,我們常用的傳統(tǒng)反應(yīng)堆動力學(xué)方程組并不具備處理上述物理現(xiàn)象的功能。這是因為,如果在系統(tǒng)中僅有少量的中子、或有隨機因素時,則在系統(tǒng)內(nèi)、各空間點上單位體積元內(nèi)中子數(shù)的多少也許不是平均值,而是一個隨機事件,系統(tǒng)中的中子數(shù)可能會非常迅速的增殖,也可能會全部死絕。
雖然,MC方法具有處理此類隨機事件的獨特優(yōu)勢,也最常見。但是,MC方法運算速度慢、人力物力成本相對較高是其主要缺點。
俄羅斯學(xué)者在這一問題研究中獨具特色,尤其在反應(yīng)堆無外源啟動的理論研究和應(yīng)用實踐中成果顯著,Ю.В.沃爾科夫1988年就提出來隨機動力學(xué)的概念。
1992年在他的《弱源條件下的反應(yīng)堆隨機動力學(xué)與核安全研究》[4]一文中簡單介紹了可將隨機微積分應(yīng)用于描述反應(yīng)堆物理中的隨機現(xiàn)象,并給出了一組點堆隨機動力學(xué)微分方程組。他指出,可以用這一方程組來描述反應(yīng)堆弱源啟動過程中中子的隨機性漲落物理現(xiàn)象、堆內(nèi)多群中子持續(xù)時間內(nèi)的隨機軌跡過程、討論中子數(shù)的隨機行為和源強對核安全影響的重要性。
雖然他沒有介紹推導(dǎo)過程,也沒有介紹求解計算方法,我們也難以直觀簡單地判斷其方程的正確性,但是作者創(chuàng)新地在點堆動力學(xué)方程中引入了隨機項的建模思路引起了我們極大的關(guān)注和興趣。
隨機微積分方法是近年來國際多學(xué)科領(lǐng)域處理各種隨機現(xiàn)象的一種熱門的數(shù)學(xué)方法,日本數(shù)學(xué)家K.Ito(伊藤清)在這一數(shù)學(xué)領(lǐng)域的成果引起了國際學(xué)術(shù)界的極大關(guān)注與贊賞。即,求解隨機微分方程最有實用意義的近似方法是伊藤清提出的Ito公式,有了Ito 公式之后,就可以計算一些基本常用的隨機微積分方程。
事實上,這一建模及數(shù)學(xué)處理方法已經(jīng)普遍在國際金融、特別是股票交易中漲落情況的預(yù)計分析、水利、工程可靠性及不確定性分析等具有隨機事件的領(lǐng)域中都看到諸多成功的研究和應(yīng)用范例[5-7]。但尚未見到有國內(nèi)學(xué)者在反應(yīng)堆物理中的應(yīng)用實踐。
我們借鑒俄羅斯學(xué)者Ю.В.沃爾科夫的思想,仔細推導(dǎo)了點堆隨機動力學(xué)微分方程組、在點堆動力學(xué)方程中引入隨機項的理論研究工作。其目的是,探索快速計算和準確評估反應(yīng)堆弱源啟動過程中的中子增殖隨機性漲落物理現(xiàn)象及不確定度,并希望引起國內(nèi)學(xué)者的興趣與關(guān)注。
本文第一節(jié),介紹了我們擬研究的物理問題,簡介了Ю.В.沃爾科夫的方程組。第二節(jié)介紹了我們從中子隨機事件的基本假定出發(fā),推導(dǎo)點堆隨機動力學(xué)微分方程組的過程和結(jié)果,并指出了Ю.В.沃爾科夫方程組的錯誤。
一般情況下,對于強中子源的核反應(yīng)堆啟動、運行、或失控時,應(yīng)用我們熟知的點堆核反應(yīng)動力學(xué)常微分方程組來描述中子的增殖過程是可信的,也不存在隨機問題,只需正確地描述并且考慮最初的無限裂變鏈出現(xiàn)的概率即可。
然而,常用的點堆反應(yīng)動力學(xué)方程在處理和描述含有弱中子源的反應(yīng)堆啟動、臨界安全分析、脈沖堆爆發(fā)脈沖實驗等物理過程時,就容易出現(xiàn)中子增殖過程中因各種隨機事件發(fā)生隨機性漲落。例如田灣核電站無源(實為燃料棒的自發(fā)裂變源)啟動的反應(yīng)性測量曲線則直觀地展現(xiàn)了弱源啟動階段反應(yīng)性與中子增殖過程中的隨機漲落現(xiàn)象[8-10],如圖1所示。
圖1 田灣核電廠無源啟動反應(yīng)性測量曲線
此時,對于這種增殖過程存在隨機性漲落,如果采用統(tǒng)計平均的處理方法是存在問題的。
再如,1960年美國學(xué)者T.F.Wimett等人公開發(fā)表了在Godiva-Ⅱ脈沖堆進行的弱源爆發(fā)脈沖實驗[3],也揭示了快堆啟動中的隨機行為和復(fù)雜的物理現(xiàn)象。其主要特征是每次爆發(fā)脈沖堆時間呈隨機性,總的爆發(fā)脈沖等待時間概率分布具有規(guī)律性,還會出現(xiàn)提前爆發(fā)脈沖現(xiàn)象,(在脈沖棒未達到預(yù)定反應(yīng)性時就已經(jīng)提前爆發(fā)脈沖)。
中國工程物理研究院物理與化學(xué)研究所也重現(xiàn)了這一現(xiàn)象[11]。
所以,在弱源條件下的核反應(yīng)堆啟動與安全運行控制時,就要特別關(guān)注反應(yīng)性、或中子數(shù)和緩發(fā)中子先驅(qū)核的隨機漲落及不確定度,因為,在弱源條件、當k=k(t)時,這種隨機漲落現(xiàn)象及不確定度將影響反應(yīng)堆的啟動時間,也關(guān)系到反應(yīng)堆的安全風(fēng)險分析的準確性。
傳統(tǒng)的點堆模型動力學(xué)方程實際上是一個數(shù)量相互作用的模型系統(tǒng),特別是中子的數(shù)量和緩發(fā)中子先驅(qū)核的數(shù)量,當物理的動力學(xué)系統(tǒng)被標示為數(shù)量過程后,在技術(shù)上可以把確定論的點堆動力學(xué)方程發(fā)展為隨機微分方程系統(tǒng),即,成為實際的隨機行為過程模型,該隨機微分方程也概括了確定論的點堆動力學(xué)方程。
如序言所述,俄羅斯學(xué)者Ю.В.沃爾科夫提出了隨機動力學(xué)的概念,介紹了可用隨機動力學(xué)微積分方程組來描述反應(yīng)堆物理中的隨機現(xiàn)象。沃爾科夫方程組形式如下:
(a)
(b)
i=2,m+1;
(c)
(d)
n=2,m+1;
(e)
(f)
n=2,m+1;
但是,由于沃爾科夫既沒有介紹推導(dǎo)過程,所以,需要我們通過推導(dǎo)和求證去判斷公式是否正確。
本節(jié)我們將推導(dǎo)和建立引入了隨機項的點堆隨機動力學(xué)微分方程組。
考慮一個十分小的時間間隔Δt,Δt時間內(nèi)有2+Q種不同事件單獨發(fā)生。因為Δt足夠小,可以假定任何兩種事件不同時發(fā)生,y1(t)表示系統(tǒng)內(nèi)中子數(shù)隨時間的變化,考慮存在兩組緩發(fā)中子β2、β3,則Δt內(nèi)有6種不同事件單獨發(fā)生,
第一種可能的貢獻,俘獲:
(1)
Δt內(nèi)俘獲中子的事件的概率為:
(2)
第二種可能的貢獻,定常中子源S在Δt內(nèi)發(fā)源射中子:
(3)
定常中子源S在Δt內(nèi)發(fā)射的概率為:
P2=SΔt
(4)
第三種可能的貢獻,Δt內(nèi)裂變釋放ν1個中子:
(5)
Δt內(nèi)裂變釋放ν1個中子的概率為:
P3=ΣfυΔty1P(ν1)
(6)
第四種可能的貢獻,Δt內(nèi)裂變釋放ν2個中子:
(7)
Δt內(nèi)裂變釋放ν2個中子的概率:
P4=ΣfυΔty1P(ν2)
(8)
第五種可能的貢獻,Δt內(nèi)先驅(qū)核y2衰變:
(9)
Δt內(nèi)先驅(qū)核y2衰變的概率為:
P5=λ2y2Δt
(10)
第六種可能的貢獻,Δt內(nèi)先驅(qū)核y3衰變:
(11)
Δt內(nèi)先驅(qū)核y3衰變的概率為:
P6=λ3y3Δt
(12)
上述式中,y1(t)表示系統(tǒng)內(nèi)中子數(shù)隨時間的變化;先驅(qū)核的份額分別為β2、β3;裂變產(chǎn)生次級中子的數(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,表達式如下:
(26)
式中,y(t)=[y1(t),y2(t),…ym+1(t)]是矢量;y1是反應(yīng)堆內(nèi)中子數(shù);
[y2(t),…ym+1(t)]是相應(yīng)的緩發(fā)中子組隊先驅(qū)核數(shù)目;
矢量ξ(t)=[ξ1(t),ξ2(t),…ξm+1(t)]由多個標準的高斯白噪聲組成;
式(26)中右側(cè)部分白噪聲的存在表示中子出生和死亡等基本過程的特征,其特征時間?L。
gij是解矩陣方程GGT=B的元素,維數(shù)為(m+1)(m+1)的矩陣,Bij的元素為:
(27)
(28)
(29)
(30)
到此,我們得到了完整的隨機中子微分方程組式(26)~方程組式(30)。
該方程組被稱為隨機點堆動力學(xué)微分方程,其調(diào)節(jié)部分與點堆動力學(xué)方程完全相同。不同的是,式中引入了隨機項∑gij(y)ξi(t),或稱為伊藤隨機項。如果矩陣方程GGT=B=0,則方程組可退化為標準的點堆動力學(xué)方程。
方程組可以描述反應(yīng)堆內(nèi)有少量的源中子和存在多群緩發(fā)中子先驅(qū)核的連續(xù)時間內(nèi)的隨機過程的軌跡。為了評估反應(yīng)堆的安全性,或發(fā)生事故時的輻射劑量,則必須弄清y(t)=[y1(t),y2(t),…ym+1(t)]的時間過程,
還必須解決中間概率下的多為擴散過程y(t)的邊界問題,該問題都可通過上述點堆隨機動力學(xué)微分方程組來描述。在引進了隨機項后,就可以來描述中子增殖過程中的隨機漲落現(xiàn)象,而不僅僅是期望值。
我們把本文得到的方程組與沃爾科夫的方程組和元素Bij進行了仔細的對比。我們發(fā)現(xiàn)不知是沃爾科夫筆誤還是印刷錯誤,在他的方程組與腳標存在以下多處錯誤:
(c)式中有bii、yi、i=1三處腳標錯誤,正確應(yīng)分別為b11、y1、i=2。
通過上述工作,我們發(fā)現(xiàn)了沃爾科夫點堆隨機動力學(xué)方程組的某些錯誤,得到了正確的方程組,這將為我們后續(xù)對方程組的求解方法研究和實驗的模擬奠定了可靠的基礎(chǔ)。我們研究認為,雖然沃爾科夫所給的方程組存在一些筆誤或錯誤,但他的物理思想是科學(xué)創(chuàng)新的。而且,隨機反應(yīng)堆動力學(xué)方程也將是今后反應(yīng)堆物理中值得研究和發(fā)展的方向之一。
本文介紹這一物理問題和方程組的推導(dǎo)過程,目的也是想引起業(yè)內(nèi)同行的興趣和關(guān)注,今后,我們將陸續(xù)給出這一方程組的求解方法研究和實驗?zāi)M的結(jié)果。
針對反應(yīng)堆物理中的隨機問題的研究,通過調(diào)研和借鑒國外學(xué)者的理念,我們嘗試著將隨機微積分理論引入于反應(yīng)堆物理,在點堆模型假定下,仔細推導(dǎo)了可描述隨機現(xiàn)象的點堆隨機動力學(xué)微分方程組,發(fā)現(xiàn)并糾正了沃爾科夫方程組的錯誤。這將為我們進行方程組的求解方法研究和實驗的模擬奠定了可靠的基礎(chǔ)。