一種快速搜索臨界棒位方法的開發與評價

李治剛，安 萍，嚴明宇,2，劉 東,3,*，蘆 韡，余紅星,2

(1.中國核動力研究設計院，四川 成都 610041；2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610041；3.中核集團核能軟件與數字化反應堆工程技術研究中心，四川 成都 610041)

目前，堆芯物理設計與計算通常采用兩步法[1]，組件計算程序生成組件庫，堆芯程序基于組件庫進行全堆芯計算。反應性控制[2-3]是堆芯核設計的重要內容，硼溶液、控制棒等是控制堆芯反應性的重要手段。以鈉冷快堆[4-5]為代表的金屬冷卻反應堆主要采用控制棒實現堆芯反應性控制，計算臨界棒位是金屬冷卻快堆物理設計中的必要內容。

在非確定論程序中，傳統的臨界搜索方法需通過多次臨界計算，具有計算量大、臨界結果不準確等缺點，一些研究機構提出了基于微擾計算的蒙特卡羅臨界搜索方法[6-7]以克服上述缺點。而在確定論程序中，目前常采用線性插值或線性外推[8]的方法獲得堆芯臨界棒位，由于控制棒的價值與控制棒插入堆芯深度并不呈簡單的線性關系[9]，通過線性插值方法往往需多次迭代計算控制棒棒位才能使堆芯達到臨界，計算效率較低。特別針對金屬冷卻快堆換料周期較長的特點，要完成整個壽期的調棒臨界計算需3～5 d，不利于堆芯設計方案的設計和優化。

為克服調棒臨界線性插值方法迭代次數多的缺點，提高臨界搜索效率，本文提出一種基于控制棒價值函數的臨界棒位快速搜索方法，并將該方法在堆芯三維中子學物理程序PBRT[10]中予以開發。對控制棒價值函數法涉及的理論及計算流程進行描述，并基于一種典型的金屬冷卻快堆堆芯布置對其進行測試，以驗證該方法的可行性和準確性。

1 理論模型

采用控制棒搜索臨界的一般計算流程如圖1所示，線性插值法和控制棒價值函數法在臨界棒位搜索環節有所差異。前者采用線性插值獲得新的臨界棒位，而后者通過內置預先得到的控制棒棒位與有效增殖因數keff的函數關系，直接計算得到新的臨界棒位。

1.1 線性插值法

假設堆芯有N組控制棒，控制棒組編號為1～N，控制棒組初始棒位H0(本文棒位均指距離堆芯底部的高度)，第i次迭代計算時的控制棒組棒位為Hi，堆芯臨界有效增殖因數為keff,critical，堆芯臨界收斂準則ε，不考慮控制棒組疊步運動。則第i次迭代后采用線性插值計算得到的新棒位Hi+1[11]為：

(1)

式中：keff,H0和keff,Hi分別為控制棒組位于初始棒位H0和Hi時的堆芯有效增殖因數。

圖1 控制棒搜索臨界計算流程Fig.1 Control rod critical computing flow

1.2 控制棒價值函數法

由于控制棒價值[12]與控制棒插入深度不呈線性關系，采用線性插值往往需多次迭代才能搜索到堆芯臨界的控制棒棒位。若能較為準確地得到控制棒價值與控制棒組棒位的函數關系，則可直接計算得到滿足臨界收斂準則的控制棒棒位?？刂瓢魞r值函數法的計算流程如圖2所示。

圖2 控制棒價值函數法計算流程Fig.2 Calculation flow of control rod value function method

Δki=keff,z(i)-keff,z(i-1)

(2)

(3)

式中：keff,z(i)和keff,z(i-1)分別為控制棒組位于z(i)和z(i-1)處堆芯的有效增殖因數；keff,bottom為控制棒組位于堆芯底部時堆芯的有效增殖因數。

假設第i組控制棒積分價值與控制棒組棒位的關系為y(i)=f(z(i))，第i組控制棒移動Δz(i)時的keff為：

keff=keff,z(i)+fi(z(i)+Δz(i))-fi(z(i))

(4)

考慮N組控制棒移動時的keff為：

keff=keff,(z(1),z(2),…,z(N))+

(5)

根據調棒順序，逐步提出控制棒組，直到堆芯keff滿足臨界收斂準則，并記錄此時各控制棒組的棒位z(i)(i=1,…,N)，即為臨界棒位。

采用最小二乘法曲線擬合[13]或切比雪夫(Chebyshev)[14-15]曲線擬合多項式擬合的方法獲取控制棒組積分價值與控制棒組棒位的函數關系fi(z(i))。

2 測試與驗證

將控制棒價值函數法應用到堆芯三維中子學程序PBRT的調棒臨界搜索模塊中，并以金屬冷卻快堆典型六邊形堆芯布置來驗證該方法在實現快速搜索臨界棒位的可行性。在堆芯中布置有3種不同價值的控制棒組，編號為1、2、3。假設共計算10個燃耗步，每個燃耗步步長為100滿功率天。

1) 控制棒組價值

采用PBRT計算得到壽期初、壽期末3組控制棒分別單獨提出堆芯的keff。圖3示出了3組控制棒微分價值和積分價值隨控制棒組棒位的變化。從圖3可看出控制棒組提出微分價值隨控制棒組的提出呈拋物線趨勢變化，在接近堆芯活性區中部位置時，引入的keff變化幅度達到最大；控制棒組提出積分價值隨控制棒組的提出呈近似S形曲線變化。提出微分價值和提出積分價值在壽期初至壽期末的變化幅度較小，因此在控制棒價值函數方法實現過程中，可假定控制棒價值在整個壽期(壽期較短)保持不變。

采用最小二乘法多項式擬合方法和切比雪夫最佳擬合方法獲得各組控制棒提出積分價值與控制棒組棒位的函數關系。圖4示出了兩種擬合方法對1號棒組的擬合結果，其他棒組的擬合規律與之基本類似。從圖4可看出，擬合階數越低，擬合函數與真實結果偏差越大；5階時，最小二乘法擬合結果優于切比雪夫法；10階時，兩種擬合方法與真實結果完全重合為S形曲線。 2) 臨界棒位搜索

將最小二乘法擬合的控制棒價值函數應用到堆芯物理設計程序的調棒臨界搜索模塊中，臨界搜索收斂準則為ε=0.001。圖5示出了控制棒價值函數法與線性插值法在第2個燃耗步時搜索臨界過程中keff的變化?？刂瓢魞r值函數法僅進行了1次迭代計算便搜索到了臨界棒位；而線性插值法需要15次迭代插值計算才滿足臨界收斂準則。

圖3 3組控制棒微分價值和積分價值隨棒位的變化Fig.3 Differential value and integral value of three sets of control rods vs. rod position

圖4 1號棒組提出積分價值函數與擬合階數的關系Fig.4 Raising integral value function of No.1 rod vs. fitting order

圖5 第2個燃耗步調棒臨界搜索過程keff變化Fig.5 keff change during process of searching critical rod position in the second burnup step

圖6示出了兩種方法計算的第2～10個燃耗步的臨界棒位和對應的keff。圖5中顯示的為2號棒組的棒位，在第2至10個燃耗步過程中，1號棒組已完全提出堆芯，3號棒組仍全部插入。從圖6可看出，兩種方法計算得到的臨界棒位的數值及趨勢基本一致，除第3個燃耗步時兩種方法的計算偏差超過4%，其他各點偏差均在2%范圍內。采用線性插值法和控制棒價值函數法計算得到的keff均在臨界搜索收斂準則范圍內。

圖6 臨界棒位和keff隨燃耗步的變化Fig.6 Critical rod position and keff vs. burnup step

3) 計算效率

圖7示出了單個燃耗步的計算時間以及燃耗步的累計計算時間。采用線性插值法的第1次調棒臨界搜索計算時間較長，后續燃耗步的計算時間基本維持在約1 h。而采用控制棒價值函數法后續燃耗步的計算時間基本維持在約0.5 h。從圖7可看出，完成10個燃耗步計算，采用控制棒價值函數法的累計計算時間較采用線性插值法的縮短了1倍，計算效率明顯提高。

3 結論

由于采用線性插值法進行調棒搜索臨界需要多次迭代計算、效率低，本文提出了一種基于控制棒價值函數的快速調棒搜索臨界方法，對該方法的理論及計算流程進行了闡述，并在堆芯物理設計程序中進行了應用和初步驗證，得出如下結論。

1) 與線性插值法相比，采用控制棒價值函數法僅進行1次迭代計算即可搜索到臨界棒位，滿足臨界搜索收斂準則，搜索效率明顯提高。

圖7 單個燃耗步的計算時間以及燃耗步的累計計算時間Fig.7 Computation time of single burnup step and cumulative computation time of all burnup steps

2) 采用控制棒價值函數法，燃耗計算時間縮短了約1倍，有助于提高堆芯物理設計的效率。

3) 通過增加計算控制棒價值的數據點和提高多項式擬合的階數，可提高控制棒價值函數法的準確度。

本文提出的控制棒價值函數法為堆芯物理設計和優化降低了計算時間，是一種合理可行的快速搜索臨界棒位方法，有效克服了傳統線性插值方法中的不足。