基于可視化方法的點堆動力學(xué)方程剛性分析

余 濤，欒秀春，劉 磊，劉少有

(哈爾濱工程大學(xué)核科學(xué)與技術(shù)學(xué)院，哈爾濱150001)

在反應(yīng)堆動態(tài)過程中，空間效應(yīng)不是主要的，如果我們需要研究中子密度隨時間的變化，此時可以應(yīng)用點堆模型[1].對應(yīng)于點堆模型下描述反應(yīng)堆中子密度及先驅(qū)核濃度變化的微分方程組稱為點堆動力學(xué)方程.由于點堆模型的簡便，目前已經(jīng)廣泛應(yīng)用于工程設(shè)計中，主要是空間效應(yīng)不太重要時，比如反應(yīng)堆、蒸汽發(fā)生器的安全事故分析.

在工程設(shè)計中需要求解點堆動力學(xué)方程，但由于物理過程中的實際情況，求解時會遇到微分方程的剛性問題，給求解帶來一定的困難，因此選擇合適的求解方法克服方程的剛性顯得尤為重要[2－5].利用可視化的方法能夠較為直觀地表現(xiàn)出動力學(xué)的剛性性質(zhì)，方便對剛性的理解.本文將利用可視化方法，形象地說明點堆動力學(xué)方程剛性問題的形成及來源.

1 點堆動力學(xué)方程

考慮緩發(fā)中子效應(yīng)后的反應(yīng)堆動態(tài)方程通常稱為點對模型動態(tài)方程，即點堆動力學(xué)方程.

不考慮中子源，緩發(fā)中子群數(shù)為6時的點堆動力學(xué)方程為:

其中:n(t)為t時刻的中子密度;Ci(t)為t時刻第i組緩發(fā)中子先驅(qū)核濃度;ρ(t)為t時刻的反應(yīng)性;Λ為中子一代時間即瞬發(fā)中子的平均壽命;λi為第i組緩發(fā)中子先驅(qū)核衰變常數(shù);βi為第i組緩發(fā)中子份額

在實際問題中，反應(yīng)性時時間的函數(shù)，如果考慮到功率和溫度對反應(yīng)性的反饋作用，ρ還是中子通量密度或n(t)的函數(shù)，因此上述方程一般來講是非線性的.本文主要探究方程的剛性問題，對方程進(jìn)行一定的簡化，變成一般的線性方程組.

2 方程的求解

上述方程的求解屬于微分方程組的初值問題，對于一階線性常系數(shù)微分方程組可以利用函數(shù)法求解，假定其解具有以下形式:

和

式(2)(3)中 ωi為特征參數(shù)，Ai、Ci，j為系數(shù).

3 方程的剛性分析

3.1 剛性系統(tǒng)

在工程中，對于一類過程的可以歸結(jié)為一般常微分方程的初值問題[6]:

1)J的所有特征值的實部Reωi＜0

則稱方程(1)為剛性方程組，r為剛性比.根據(jù)Shempine和Gear的觀點，如果矩陣J的所有特征值ωi的實部都不是很大的正數(shù)，同時至少有一個特征值是很大的負(fù)數(shù)，則該系統(tǒng)也視為具有剛性的.

3.2 點堆方程的剛性

考慮在熱堆中輸入正階躍反應(yīng)性的過程，Λ=0.000 5 s，βi:0.000 285、0.001 598、0.00 141、0.003 052、0.000 96、0.000 195;λi:0.012 7、0.031 7、0.105、0.311、1.4、3.87;輸入的正階躍反應(yīng)性:ρ=0.001，計算時間為 100 s.

根據(jù)圖解法可求得:ωi:0.018 2、－0.013 6、－0.059 8、－0.183、－1.005、－2.875、－55.6，設(shè)定ω1＜ ω2＜ …… ＜ ω7，f(i)=Aieωit(i=1，2，…，7).在t=0 s時，反應(yīng)堆中輸入階躍反應(yīng)性，fi隨時間增加，圖1中f(1)從t=0 s時刻反應(yīng)性輸入后持續(xù)增加，以指數(shù)形式增大，變化趨勢較為明顯;f2、f3、f4、f5、f6等 5 條曲線變化趨勢幾乎一致，在圖中部分曲線近乎重合，同時圖中基本看不出這幾條曲線的變化趨勢，近乎為一條直線;f7曲線在t=0 s時刻后有一個增大過程，但圖1中無法分辨變化幅度與變化變化時間.圖2為f1～f7在不同坐標(biāo)系下的變化曲線，圖中明顯可以看出f1的值比其他幾項大，所以造成了圖1中f2～f7的幅值變化無法清楚地顯示.

圖1 反應(yīng)性正階躍輸入時0～100 s內(nèi)中子密度變化

圖2 f1～f7在不同坐標(biāo)系下的變化曲線

圖3、4 分別為f2、f3、f4和f5、f6的變化.對比圖3、4中各項與圖1中各對應(yīng)項，圖3較為清楚地展現(xiàn)了各項的變化規(guī)律，f2、f3、f4項在100 s內(nèi)持續(xù)明顯的增大.在圖2中，延長計算時間后發(fā)現(xiàn)，f1持續(xù)增大，f2在t=450 s后趨于0，f3在t=100 s后趨于0，f4在t=30 s后趨于0;圖4中可以看出f5、f6在反應(yīng)性輸入后一段時間內(nèi)增大，然后趨于0，趨于0的時間點分別為:5、2 s，而在圖1中則無法看出.

圖5為f7在0～1s內(nèi)變化曲線，在0～0.5 s內(nèi)，該項由－0.179變化到0左右，0.1 s之后則基本為0，變化趨勢與f2～f6相同.由于該項的增長時間較短，因此在圖1中無法分辨出增長時間段長度.從上述計算結(jié)果圖示可以看出，式(2)中右邊各項在反應(yīng)性輸入后均隨時間增加，第1項持續(xù)明顯增加，沒有趨于一個穩(wěn)定值，其余6項也隨時間增加，但在結(jié)果一段時間后趨于穩(wěn)定值零，各項趨于穩(wěn)定值的時間依次減小.對于圖1中的曲線，由于各曲線對時間的響應(yīng)速度不一致，圖中無法清楚顯示各項的變化，因此需要在不同圖示中顯示各項的變化，這正是方程組剛性性質(zhì)存在所造成的.

圖5 f7的變化曲線

對于正階躍反應(yīng)性輸入時，式(2)中:A1＞0，ω1＞0;Ai＜0，ωi＜0(i=2，…7)，從而造成上述f1變化趨勢與其他各項不同的現(xiàn)象，在f2～f7曲線中，由于各Ai、ωi之間數(shù)量級的差異較大使得圖1中顯示的曲線近乎重合.圖5中，f7的變化時間為0 ～0.1s，在 0.1s內(nèi)，f7有一個較大的增大過程，在相同時間段內(nèi)，f1～f6的變化量較小，說明在0～0.1s時間段內(nèi)，中子密度的變化主要是由f7貢獻(xiàn)的，而在0.1s后，f7趨于零.圖1中顯示的則是在0s時刻的一個階躍，即曲線斜率趨于+∞，說明在一個比較大的時間尺度下，f7無法正常體現(xiàn)出來.對應(yīng)于在式(1)的數(shù)值求解過程中，如果選取的時間步長過大，則對于方程的穩(wěn)定性無法保證，因此由于f7的存在需要在數(shù)值求解時選取較小的時間步長，保證計算方法的穩(wěn)定性以及計算結(jié)果的正確性.但是如果選擇較小的計算步長，則對計算效率和計算時間產(chǎn)生影響.由于f7對中子密度時間響應(yīng)的貢獻(xiàn)只存在于0.1 s內(nèi)，0.1 s后基本不需要考慮f7的影響，也就不需要采用因該項而決定的時間步長，從而造成計算效率的浪費.過小的計算時間步長在進(jìn)行較長計算時間的計算時意味著較長的工時.在工程設(shè)計中，上述選擇計算時間步長限制的矛盾正是系統(tǒng)剛性性質(zhì)的表現(xiàn).

由于剛性微分方程的求解特殊性，因此在工程設(shè)計中盡量避免剛性微分方程，但由于實際過程的限制，實際過程中的多數(shù)方程都具有一定的剛性，其中一些的剛性還較大.對于式(1)，根據(jù)剛性系統(tǒng)的定義，在 ρ=0.001階躍輸入時，方程符合Shempine和Gear對于剛性系統(tǒng)的觀點;在ρ=－0.001階躍輸入時，方程剛性比r為1 000，符合剛性系統(tǒng)的定義.在剛性系統(tǒng)的定義中，判斷是否是剛性系統(tǒng)取決于特征值 ω，對于點堆方程由下式[7]:

ω1＞ －λ1＞ω2＞… ＞ω6＞ －λ6＞ω7＞l－1(5)其中:l為中子壽命，根據(jù)式(5)可知，由于反應(yīng)堆緩發(fā)中子先驅(qū)核衰變常數(shù)的限制使得方程(1)的各個特征值ωi之間相差較大，對于熱堆，λ1大約為10－2數(shù)量級，λ6為 101數(shù)量級，l－1為 103數(shù)量級，因此在負(fù)階躍反應(yīng)性輸入時，ωi＜0，此時方程剛性比為102～104數(shù)量級，屬于剛性系統(tǒng).

4 結(jié)語

在可視化方法的基礎(chǔ)上，分析了點堆動力學(xué)方程的剛性性質(zhì)，方程剛性是由核反應(yīng)堆的物理性質(zhì)決定的，即堆內(nèi)緩發(fā)中子先驅(qū)核的衰變常數(shù)以及中子一代時間之間相差數(shù)個數(shù)量級，從而使方程具有較大的剛性.

[1]謝仲生.核反應(yīng)堆物理分析[M].西安:西安交通大學(xué)出版社，2004.

[2]陳 東.求解點堆中子動力學(xué)方程的三角插值法[J].原子能科學(xué)技術(shù)，2010，44(10):1195 －1200.

[3]胡大璞.克服點堆中子動力學(xué)方程剛性的新方法——端點浮動法[J].核動力工程，1993，14(2):122－128.

[4]黎浩峰.用高階泰勒多項式積分方法求解點堆中子動力學(xué)方程[J].原子能科學(xué)技術(shù)，2008，42(增刊):162－168.

[5]NAHLA A A.Taylor’s series method for solving the nonlinear point kinetics equations[J].Nuclear Engineering and Design，2011，241(5):1592－1595.

[6]黃祖洽.核反應(yīng)堆動力學(xué)基礎(chǔ)[M].北京:北京大學(xué)出版社，2007.

[7]謝仲生.核反應(yīng)堆物理數(shù)值計算[M].北京:原子能出版社，1997.