時空中子動力學計算的剛性限制法

魏明哲，李云召，何明濤

(西安交通大學 核科學與技術學院，陜西 西安 710049)

時空中子動力學計算的剛性限制法

魏明哲，李云召*，何明濤

(西安交通大學 核科學與技術學院，陜西 西安 710049)

隨著中子動力學計算的精度要求以及計算條件的不斷提高，如何更好地處理三維時空中子動力學方程的剛性問題對于提高計算效率具有非常重要的工程應用價值。本文應用剛性限制法處理方程組的時間變量，并與空間變量處理方法變分節塊法相耦合建立了相應的時空中子動力學計算模型。基準題的數值計算結果表明，剛性限制法可與包括變分節塊法在內的各種穩態計算方法耦合以求解時空多群中子動力學問題，且取得了很好的計算精度，允許較大的計算時間步長，以及擁有很好的計算效率。

剛性限制法；點堆中子動力學；時空中子動力學；變分節塊法

剛性限制法(SCM)是一種中子動力學的時間變量處理方法，最早由美國西屋公司的Chao等[1-2]于20世紀80年代提出。其主要思想是將先驅核濃度控制方程的剛性消除掉，使剛性僅限制在中子密度方程中，而該方程可用解析法求解，從而使整個系統的求解時間步長增大。

剛性限制法的實際應用價值在于快速求解時空中子動力學模型，對于該模型，剛性限制法通過引入動力學頻率，將三維多群時空動力學問題轉換為三維多群穩態問題，因而緩解動力學問題的“剛性”，使與剛性限制法耦合的穩態程序可在大時間步長下快速計算瞬態變化，并擁有足夠的穩定性、準確性和計算效率。

美國西屋公司的趙榮安博士于1984年首次將剛性限制法應用于點堆動力學問題的求解中[1]；并于1989年，在西屋公司的壓水堆節塊程序SUPER-NOVA(SPNOVA)的基礎上，開發出了三維動力學計算程序SPNOVA-K[2]。之后，很多學者對此方法進行了研究與發展：2007年日本三菱公司的AOKI等[3]用剛性限制法耦合三維擴散程序ANC得到其動力學版本ANCK，并用ANCK與MIDAC耦合來準確模擬反饋效應；2012年上海核工程研究設計院的司勝義開發了通用的剛性限制法計算程序UASCM[4]，可直接與多種穩態擴散或輸運程序耦合應用，計算動力學問題。本文基于剛性限制法的主要思想建立計算模型，應用變分節塊法(VNM)計算時空中子動力學問題，并對計算效率進行相應優化。

1 剛性限制法的理論模型

瞬發中子動力學頻率：

(1)

緩發中子動力學頻率：

(2)

其中：n(t) 為時刻t的中子密度，cm-3；Ci(t) 為時刻t的第i組緩發中子先驅核濃度，cm-3。

1.1 點堆中子動力學方程的求解模型

點堆中子動力學方程的一般形式[1]為：

(3)

(4)

其中：Λ為瞬發中子的平均壽命；ρ為當前時刻反應性；λi為第i組緩發中子先驅核衰變常量；βi為第i組緩發中子的份額，β=∑βi。

將式(1)代入式(3)，可得：

(5)

將上式代入式(4)，并消除中子密度項，可得：

(6)

采用全隱式向后差分方法求解式(6)，得到：

(7)

(8)

其中：

(9)

(10)

將求得的緩發中子先驅核濃度代入式(3)，得：

(11)

假設在每個時間步長內，緩發中子先驅核濃度不變，用解析法求解該方程：

(12)

將得到的n(tn+1)與Ci(tn+1)代入式(3)計算dn(t)/dt，并更新動力學頻率，得：

(13)

重復以上過程直至動力學頻率收斂，該時間步計算結束。

1.2 時空中子動力學的求解模型

在笛卡爾坐標系中，分群形式的三維時空擴散中子動力學方程組[1]可描述為：

(14)

(15)

其中：vg為第g群中子的(群)平均速度，cm/s；r=(x,y,z)為空間位置；φg(r,t)為在t時刻r位置處第g群的中子通量密度，cm-2·s-1；Dg(r,t)為第g群的擴散系數；Σt,g(r,t)為第g群中子的總宏觀截面，cm-1；χig為第i組緩發中子出現在g群中的概率；χg為裂變中子進入第g群的份額；Σg′-g(r,t)為t時刻r位置處g′群向g群的轉移宏觀截面，cm-1。

(16)

引入動力學頻率，將式(14)和(15)聯立，并消去中子通量密度方程中的緩發中子先驅核濃度項：

(17)

將上式化為穩態形式的三維多群中子擴散方程：

χg(1-β)Q′(r,t)+Σg′-g(r,t)φg′(r,t)

(18)

其中：

(19)

(20)

(21)

(22)

利用穩態多群中子擴散計算程序即可求解上述方程組。

(23)

將其代入式(15)，可求得：

(24)

得到系數A后，便可求得中子通量密度的絕對值：

(25)

按照定義更新動力學頻率的值：

(26)

類似地，重復以上過程直至動力學頻率收斂，該時間步計算結束。

時空剛性限制法通過引入動力學頻率，將三維多群時空動力學問題轉換成三維多群穩態擴散問題，因而緩解動力學問題的“剛性”，使與剛性限制法耦合的穩態程序可在大時間步長下快速計算瞬態變化。

2 數值結果及其分析

本文應用上述理論模型首先編制了剛性限制法點堆計算程序SCM-Point，并將其嵌入改進準靜態計算的點堆方程求解過程中。其次在變分節塊法程序Violet的基礎上實現其時空動力學版本Violet-K。采用TWIGL Seed Blanket瞬態基準題對上述兩個部分分別進行驗證，并對計算結果進行分析研究。

TWIGL Seed Blanket瞬態基準題是一個二維兩群的瞬態問題。堆芯是一邊長為160 cm的正方形，由兩種燃料材料組成并分為3個區域，無反射層和反饋。其中，1號區域與2號區域材料相同，1號區域模擬了控制棒插入的情況。問題詳細描述參見文獻[4]。

該基準題由兩組問題組成，分別是燃料區1的熱群宏觀吸收截面發生階躍變化和線性變化引起的瞬態過程，初始時刻堆芯相對功率為1.0，瞬態過程計算到0.5 s。

2.1 點堆動力學計算數值結果

選取改進準靜態計算過程中0～10 ms的時間段。通量密度以對于初始通量密度的比值即相對通量密度表示。在線性變化和階躍變化過程，反應性引入量分別為18.48 pcm和383.06 pcm。在此時間段內分別取時間步數1、5、10、50、100、1 000計算時間節點10 ms處的通量密度，采用SCM-Point和隱式差分方法得到的數值計算結果列于表1。

通過與傳統的隱式差分方法進行比較后發現，欲達到相同計算精度，剛性限制法可以采用較大的時間步長，本例中剛性限制法的步長可以是隱式差分方法步長的數十倍。數值驗證表明，點堆模型中，剛性限制法能夠更好地處理方程的剛性，使得在相同的計算精度下，剛性限制法較隱式差分方法收斂更快，可以在保持穩定性的前提下，增大計算時間步長，從而提高計算效率。

表1 TWIGL Seed Blanket瞬態基準問題10 ms時刻相對通量密度Table 1 Relative flux density of TWIGL Seed Blanket benchmark problem at 10 ms

2.2 時空動力學計算數值結果

使用Violet-K程序選取1/4堆芯進行計算，節塊劃分為10×10，即x軸方向節塊數為10，y軸方向節塊數為10。在20 ms和10 ms兩種時間步長下計算堆芯相對功率隨時間的變化，并列出計算時間，參考解取自文獻[5]，計算結果列于表2。計算過程中，動力學頻率的更新值取前一次迭代與后一次迭代的中間值，外迭代條件取10-6，有效增殖因數收斂條件取10-6，由此得到單步的迭代計算次數約為3次。

表2 基于Violet-K程序下TWIGL Seed Blanket瞬態基準問題堆芯相對功率隨時間的變化Table 2 Variation of core relative power with time in TWIGL Seed Blanket benchmark problem based on Violet-K code

表3列出基于3種穩態程序Violet-K、MGSNM與MGNEM的剛性限制法對于TWIGL Seed Blanket瞬態基準問題的數值結果比較，后兩個同類程序數據來源于文獻[6]。在相同網格劃分條件下，分別對20 ms和10 ms的時間步長，比較3種程序下堆芯相對功率隨時間的變化。

數值結果表明，在20 ms和10 ms兩種時間步長下，堆芯相對功率的計算結果與參考值均符合得較好，尤其是當步長取20 ms時，仍可保證很好的計算精度，而隱式差分方法對于此問題為保證計算精度通常選取的時間步長為1 ms或2.5 ms。可見時空剛性限制法可較大地增大計算步長，在提高計算效率的同時保證準確率和穩定性。

通過比較基于3種穩態程序下剛性限制法關于堆芯相對功率隨時間的變化計算的結果發現，在同類的時空動力學計算程序中，本文所編制的Violet-K程序具有較良好的精度和穩定性，計算結果準確可靠。

表3 基于3種穩態程序下TWIGL Seed Blanket瞬態基準問題堆芯相對功率隨時間的變化Table 3 Variation of core relative power with time in TWIGL Seed Blanket benchmark problem based on different codes

3 結論

數值結果驗證表明，剛性限制法能更好地處理方程的剛性問題，使得在相同的計算精度下，其可在保持穩定性的前提下，增大計算時間步長，從而提高計算效率。剛性限制法可與包括變分節塊法的各種穩態計算方法耦合以求解時空多群中子動力學問題，且取得了很好的計算精度，其允許較大的計算時間步長，擁有很好的計算效率。

[1] CHAO Y, ATTARD A. A resolution of the stiffness problem of reator kinetics[J]. Nuclear Science and Engineering, 1985, 90: 40-46.

[2] CHAO Y, HUANG P. Theory and performance of the fast-running multidimensional pressurized water-reactor kinetics code, SPNOVA-K[J]. Nuclear Science and Engineering, 1989, 103: 415-419.

[3] AOKI S, SUEMURA T, OGAWA J, et al. The verification of 3 dimensional nodal kinetics code ANCK using transient benchmark problems[J]. Journal of Nuclear Science and Technology, 2007, 44(6): 862-868.

[4] SI Shengyi. Algorithm development and verification of UASCM for multi-dimension and multi-group neutron kinetics model[C]∥PHYSOR 2012: Adcances in Reactor Physics-Linking Research, Industry, and Education. USA: [s. n.], 2012.

[5] BANDINI B R. A three-dimensional transient neutronics routine for the TRAC-PF1 reactor thermal hydraulic computer code[D]. Pennsylvania: Pennsylvania State University, 1990.

[6] 吳宏春. 中子擴散方程及其數值方法[M]. 西安:西安交通大學出版社,2005.

Stiffness Confinement Method for Space-time Neutron Kinetics Calculation

WEI Ming-zhe, LI Yun-zhao*, HE Ming-tao

(SchoolofNuclearScienceandTechnology,Xi’anJiaotongUniversity,Xi’an710049,China)

With the improvement of computational ability and the demand for accuracy requirement of neutron kinetics, it has great value in engineering application to deal with the stiffness of three-dimensional space-time neutron diffusion kinetics better. Based on the stiffness confinement method (SCM) and the variational nodal method (VNM), a new combination of neutron kinetics model was developed. It is demonstrated by the numerical results that the SCM can handle the stiffness problem better with larger time step and higher computational efficiency, and be combined with the VNM.

stiffness confinement method; point neutron kinetics; space-time neutron kinetics; variational nodal method

2014-06-29;

2014-08-25

國家自然科學基金資助項目(91226106)；長江學者和創新團隊發展計劃資助項目(IRT1280)；上海核工程研究設計院資助項目

魏明哲(1992—)，男，遼寧鐵嶺人，核工程與核技術專業

*通信作者：李云召，E-mail: yunzhao@mail.xjtu.edu.cn

TL329.2

A

1000-6931(2015)10-1828-05

10.7538/yzk.2015.49.10.1828