蔡光明，阮良成

（中國核電 福建福清核電有限公司，福建 福清 350318）

點(diǎn)堆中子動(dòng)力學(xué)方程描述反應(yīng)堆中子密度和反應(yīng)性之間隨時(shí)間的變化。求解中子密度隨時(shí)間的變化，對(duì)于反應(yīng)堆的事故分析、仿真、運(yùn)行及物理試驗(yàn)都有重要的作用。

幾十年來，點(diǎn)堆中子動(dòng)力學(xué)方程求解一直是個(gè)研究熱點(diǎn)。由于點(diǎn)堆中子動(dòng)力學(xué)方程是個(gè)剛性方程[1－3］，即方程中的幾個(gè)時(shí)間常數(shù)之間相差幾個(gè)數(shù)量級(jí)，因此準(zhǔn)確、快速、穩(wěn)定地求解方程是困難的。得益于現(xiàn)代計(jì)算機(jī)技術(shù)的進(jìn)步，本文直接采用代中子時(shí)間計(jì)算法求解點(diǎn)堆中子動(dòng)力學(xué)方程，經(jīng)過基準(zhǔn)例題和動(dòng)態(tài)—逆動(dòng)態(tài)對(duì)比計(jì)算，驗(yàn)證了模型、程序計(jì)算的準(zhǔn)確性和穩(wěn)定性，而計(jì)算時(shí)間也是可接受的。

1 代中子時(shí)間計(jì)算法

當(dāng)反應(yīng)堆處于動(dòng)態(tài)時(shí)，為了簡化模型，假設(shè)堆內(nèi)中子注量率各處變化與空間位置無關(guān)，即可得到假定“點(diǎn)堆，單群，6組緩發(fā)中子”的點(diǎn)堆中子動(dòng)力學(xué)方程：

其中：n（t）為與時(shí)間相關(guān)的中子密度；βi為第i組緩發(fā)中子份額；I為緩發(fā)中子平均價(jià)值因子；βieff為第i組有效緩發(fā)中子份額；βeff為有效緩發(fā)中子份額；λi為第i組緩發(fā)中子先驅(qū)核的衰變常數(shù)；ci（t）為第i組先驅(qū)核密度；l0為瞬發(fā)中子平均壽命；l為瞬發(fā)中子平均代時(shí)間；S為中子源；keff是反應(yīng)堆有效增殖系數(shù)。

由于計(jì)算機(jī)技術(shù)的發(fā)展，本計(jì)算模型的中心思想是直接采用瞬發(fā)中子平均代時(shí)間作為時(shí)間增量來計(jì)算中子密度的變化，以計(jì)算時(shí)間為代價(jià)來解決方程的剛性問題，來保證點(diǎn)堆中子動(dòng)力學(xué)方程解的準(zhǔn)確性和穩(wěn)定性。因此，令dt＝l，則方程（1）、（2）分別轉(zhuǎn)化為：

則下一代的中子密度、第i組先驅(qū)核密度分別為：

其中，反應(yīng)性ρ定義為

2 程序與驗(yàn)證

根據(jù)第1節(jié)所述的模型，編制了計(jì)算機(jī)程序 KNGT（Kinetic Code of Neutron Generation Time method），并用參考文獻(xiàn)中的基準(zhǔn)例題進(jìn)行了計(jì)算比對(duì)。以下是關(guān)于程序的介紹與例題比對(duì)結(jié)果。

2.1 程序編制說明

代中子時(shí)間計(jì)算法求解點(diǎn)堆中子動(dòng)力學(xué)方程計(jì)算機(jī)程序KNGT，采用C＋＋語言設(shè)計(jì)，編譯的目標(biāo)代碼執(zhí)行效率高，可以充分發(fā)揮計(jì)算機(jī)的速度。另外所有的浮點(diǎn)變量全部定義為浮點(diǎn)雙精度，盡量避免賦值誤差和計(jì)算產(chǎn)生的截?cái)嗾`差[3］。

2.2 基準(zhǔn)例題驗(yàn)證

例1，計(jì)算熱堆輸入階躍正反應(yīng)性后，堆內(nèi)中子密度n（t）隨時(shí)間t的變化。其中子動(dòng)力學(xué)參數(shù)為：βieff＝0.000 266、0.001 491、0.001 316、0.002 849、0.000 896、0.000 182，λi＝0.012 7、0.031 7、0.115、0.311、1.4、3.87s－1，i＝1～6，l＝2.0E－5s；其他參數(shù)為：ρ0＝0，ρ＝0.003，S＝0s－1。計(jì)算結(jié)果列于表1。

表1 階躍正反應(yīng)性輸入后的熱堆中子密度變化Table1 Neutron density variation in thermal reactor after the input of positive step reactivity

續(xù)表

例2，計(jì)算熱堆輸入大階躍正反應(yīng)性ρ＝0.007后，堆內(nèi)中子密度n（t）隨時(shí)間t的變化，其他參數(shù)同例1。計(jì)算結(jié)果列于表2。

例3，計(jì)算熱堆輸入階躍負(fù)反應(yīng)性ρ＝－0.007后，堆內(nèi)中子密度n（t）隨時(shí)間t的變化，其他參數(shù)同例1。計(jì)算結(jié)果列于表3。

表2 大階躍正反應(yīng)性輸入后的熱堆中子密度變化Table2 Neutron density variation in thermal reactor after the input of great positive step reactivity

表3 階躍負(fù)反應(yīng)性輸入后的熱堆中子密度變化Table3 Neutron density variation in thermal reactor after the input of negative step reactivity

例4，計(jì)算熱堆輸入線性反應(yīng)性后，堆內(nèi)中子密度n（t）隨時(shí)間t的變化，計(jì)算結(jié)果列于表4。線性反應(yīng)性變化率為0.000 7s－1，其他參數(shù)同例1。

例5，計(jì)算堆內(nèi)引入正弦反應(yīng)性后，堆內(nèi)中于密度n（t）隨時(shí)間t的變化。ρ＝ρ0sin（πt／T0），其中，ρ0和T0為系數(shù)。計(jì)算結(jié)果列于表5。βeff＝0.007 9，λ＝0.077s－1，l＝10－8s。

表4 線性反應(yīng)性輸入后的熱堆中子密度變化Table4 Neutron density variation in thermal reactor after the input of liner reactivity

表5 復(fù)雜反應(yīng)性輸入后的反應(yīng)堆中子密度峰值Table5 Neutron density peak in thermal reactor after the input of complex reactivity

從以上例1～3幾種計(jì)算方法的比較，本文方法的計(jì)算結(jié)果與精確解（埃爾米特法，h＝0.000 1s[4］）比其他方法更吻合；且在計(jì)算時(shí)間上也并無劣勢。例4中本文方法的計(jì)算結(jié)果也比其他方法更吻合。例5反映了本文方法對(duì)復(fù)雜反應(yīng)性計(jì)算的穩(wěn)定性。

3 其他驗(yàn)證

上一節(jié)主要用基準(zhǔn)例題作了計(jì)算比較，也是各位學(xué)者一般采用的方法。本節(jié)將用動(dòng)態(tài)－逆動(dòng)態(tài)的計(jì)算結(jié)果互相比較，來做一步驗(yàn)證。

3.1 動(dòng)態(tài)－逆動(dòng)態(tài)計(jì)算比較（反應(yīng)性—中子密度—反應(yīng)性）

例6，表6是模擬動(dòng)態(tài)刻棒法DRWM的反應(yīng)性ρdrwm變化，所用的動(dòng)力學(xué)參數(shù)同例1，ρ0＝－0.01×10－5，S＝5.0E－10s－1。n是隨時(shí)間進(jìn)行的點(diǎn)堆動(dòng)力學(xué)動(dòng)態(tài)計(jì)算結(jié)果；根據(jù)n用ODRM[5］程序進(jìn)行逆動(dòng)態(tài)計(jì)算得到反應(yīng)性ρodrm，n的計(jì)算時(shí)間間隔為0.1s。

ρodrm－ρdrwm是反應(yīng)性通過動(dòng)態(tài)計(jì)算得到中子密度再通過逆動(dòng)態(tài)計(jì)算得到的反應(yīng)性的比較，這個(gè)偏差最大值為2.4×10－5，相對(duì)偏差小于0.3%；±50×10－5范圍內(nèi)，偏差最大值為0.5×10－5，相對(duì)偏差最大1%，反應(yīng)性吻合很好。由此分析可以得出KNGT的計(jì)算是正確的，準(zhǔn)確的；同時(shí)另一方面也驗(yàn)證了逆動(dòng)態(tài)程序ODRM模型與程序的正確性。

表6 動(dòng)態(tài)－逆動(dòng)態(tài)計(jì)算反應(yīng)性比較Table6 Comparison between kinetic－inverse kinetic contrast calculation on reactivity

續(xù)表

3.2 逆動(dòng)態(tài)－動(dòng)態(tài)計(jì)算比較（中子密度—反應(yīng)性—中子密度）

例7，圖1中I是某反應(yīng)堆稀釋刻棒的中子密度變化實(shí)測曲線，ρ是用ODRM程序進(jìn)行逆動(dòng)態(tài)計(jì)算得到反應(yīng)性，n是用KNGT進(jìn)行的點(diǎn)堆動(dòng)力學(xué)動(dòng)態(tài)計(jì)算得到的。I、ρ、n計(jì)算的時(shí)間間隔為1s。

n－I是中子密度通過逆動(dòng)態(tài)計(jì)算再通過動(dòng)態(tài)計(jì)算的比較，這個(gè)偏差最大值為0.082，最大相對(duì)偏差小于0.11%；中子密度吻合很好，圖1中子密度n、I曲線重疊無法區(qū)分。由此分析可以再次得出KNGT的計(jì)算是正確的、準(zhǔn)確的；同時(shí)另一方面也再次驗(yàn)證了逆動(dòng)態(tài)程序ODRM模型與程序的正確性。

圖1 刻棒中的中子密度與反應(yīng)性的變化Fig.1 Variation of neutron density and reactivity in rod measurement test

4 小結(jié)

對(duì)于用代中子時(shí)間法編制的程序求解點(diǎn)堆中子動(dòng)力學(xué)方程，經(jīng)過例題驗(yàn)證及動(dòng)態(tài)－逆動(dòng)態(tài)的驗(yàn)證，有如下結(jié)論：

（1）經(jīng)過基準(zhǔn)例題比較計(jì)算驗(yàn)證，模型是正確的，程序KNGT是正確的，計(jì)算結(jié)果是準(zhǔn)確的。

（2）由于模型使用了最小的時(shí)間參數(shù)，未采用特別的假設(shè)，因此模型的計(jì)算結(jié)果是穩(wěn)定的，可用于各種復(fù)雜動(dòng)態(tài)計(jì)算。

（3）在當(dāng)前計(jì)算機(jī)硬件條件下，程序運(yùn)行是足夠快的，可用于熱堆的實(shí)時(shí)仿真計(jì)算。

（4）使用動(dòng)態(tài)－逆動(dòng)態(tài)相互計(jì)算比較，更進(jìn)一步驗(yàn)證了動(dòng)態(tài)程序KNGT的計(jì)算正確性與穩(wěn)定性；另一方面也驗(yàn)證了逆動(dòng)態(tài)程序ODRM計(jì)算的正確性。

需要另外指出的是，例6中中子源S≠0，這在參考文獻(xiàn)中鮮有類似的基準(zhǔn)例題驗(yàn)證。除了本文的方法驗(yàn)證之外，希望得到其他學(xué)者的驗(yàn)證。

[1]廖忠岳，李鐵軍，劉長欣 .反應(yīng)堆中子動(dòng)力學(xué)方程仿真算法及其應(yīng)用[J］.核動(dòng)力工程，1995，16（4）：306－311.

[2]黎浩峰，等 .用高階泰勒多項(xiàng)式積分方法求解點(diǎn)堆中子動(dòng)力學(xué)方程[J］.原子能科學(xué)技術(shù)，2008，42：162－168.

[3]經(jīng)滎清，等 .點(diǎn)動(dòng)態(tài)學(xué)方程幾種解法的比較[J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），1995，35（3）：7－12.

[4]陳東，劉振華 .求解點(diǎn)堆中子動(dòng)力學(xué)方程的三角插值法[J］.原子能科學(xué)技術(shù)，2010，44（10）：1195－1120.

[5]蔡光明 .落棒法測量所有控制棒全插時(shí)的棒價(jià)值[J］.核科學(xué)與工程，2005，25（2）：154－158.