鈉冷快堆系統分析程序FR-Sdaso衰變熱計算模型開發與驗證

楊 軍，賈鴻玉，楊曉燕，王曉坤

(中國原子能科學研究院 反應堆工程技術研究部，北京 102413)

反應堆長時間運行后停堆，裂變功率迅速下降，而堆芯衰變熱則緩慢下降，堆芯余熱在停堆初期由裂變功率和衰變熱構成，停堆后期主要由衰變熱構成。停堆余熱根據產生方式主要分為以下幾項：1) 裂變產物及裂變產物中子活化后的衰變熱；2) 放射性錒系核素的衰變熱；3) 緩發中子和自發裂變導致的殘余裂變熱；4) 結構材料中子活化后的衰變熱[1]。雖然衰變熱一般只占到反應堆額定功率的百分之幾，但若不及時導出，將導致燃料元件因過熱而燒毀，因此反應堆停閉時的余熱導出是反應堆設計時需考慮的重要安全功能之一[2]。

系統分析程序是反應堆工況設計分析和事故分析的重要計算工具。目前國際上發展鈉冷快堆的國家大都開發了快堆系統分析程序，如美國阿貢國家實驗室(ANL)開發的SAS4A/SASSYS-1程序[3]、法國原子能委員會(CEA)開發的OASIS程序[4]、韓國原子能研究院(KAERI)開發的MARS-LMR程序[5]、西安交通大學開發的THACS程序[6]、華北電力大學開發的SAC-CFR程序[7]、哈爾濱工程大學開發的THPCS程序[8]等。

FR-Sdaso程序為中國原子能科學研究院針對大型池式鈉冷快堆開發的系統瞬態分析程序，其計算范圍涵蓋鈉冷快堆主熱傳輸系統，可分析鈉冷快堆電廠在正常運行工況及各類事故工況下的瞬態響應[9]，已用于示范快堆工況設計及安全分析。本文開發FR-Sdaso程序中所使用的衰變熱計算模型，并對模型進行驗證。

1 衰變熱計算方法

停堆后衰變熱的計算一般有兩種方法[10-11]：累計方法和量熱法。

累計方法單獨處理堆內數百種裂變產物中的每一種裂變產物的衰變熱，然后相加求得反應堆總的衰變熱。該方法主要依賴于裂變產物數據的正確性，需要專用的計算軟件，典型如美國橡樹嶺國家實驗室編制的ORIGEN程序。累計方法的上述特點，導致其難以應用于系統分析程序。

量熱法使用停堆后的衰變熱積分實驗曲線，進行指數多項式擬合，然后利用擬合公式計算。典型的量熱法計算公式[12]如下：

Pd(τ,T)=4.1×1011P[τ-0.2-(τ+T)-0.2]

(1)

其中：Pd(τ,T)為停堆后的衰變熱，MeV/s；P為停堆前的功率，W；τ為停堆后的時刻，s；T為停堆前維持功率P連續運行的時間，s。

一方面量熱法只能單獨計算衰變熱，不能考慮瞬態過程中裂變功率對衰變熱的影響，另一方面目前尚無成熟的針對鈉冷快堆的量熱法計算公式。因此，本文開發一種可用于鈉冷快堆系統分析程序模擬衰變熱的集總計算方法，可考慮裂變功率和功率歷史對衰變熱的影響。

1.1 衰變熱計算模型

某種裂變產物或錒系核素的濃度可近似計算[12]如下：

(2)

其中：hj(t)為t時刻第j種裂變產物或錒系核素的濃度，m-3；γj為第j種裂變產物或錒系核素的產額；Σf為總裂變截面，m-1；φ(t)為t時刻中子通量密度，m-2·s-1；λj為第j種裂變產物或錒系核素的衰變常量，s-1。實際應用中可用式(2)近似模擬1組裂變產物或錒系核素濃度的變化。

在式(2)兩邊同乘以第j種裂變產物或錒系核素每次衰變釋放出的能量Edj(J)，有：

(3)

令Edjhj(t)=Hj(t)，有：

(4)

對式(4)做如下變換：

(5)

其中：Ef為單次裂變釋放的能量，J；EfΣfφ(t)=N(t)為裂變功率密度，W/m3。

令Edjγj/Ef=Ej為無量綱數，則式(5)變換為：

(6)

其中：裂變功率密度N(t)可由點堆模型計算；Ej為常數，表示各核素對應的功率份額；λjHj(t)為第j種裂變產物或錒系核素t時刻的衰變熱密度。應用式(6)即可計算衰變熱[13]。顯然，瞬態過程中裂變功率將對衰變熱產生影響。

1.2 初始條件確定

反應堆停堆后的衰變熱受反應堆運行功率歷史的影響，因此初始化過程需考慮不同功率歷史對衰變熱計算的影響。

假設某反應堆的功率歷史可簡化為圖1所示的N段，每個分段內反應堆總功率為常數，分別為P1，P2，P3，…，PN，每個時間段的時刻分割點分別為T1，T2，T3，…，TN。

圖1 反應堆功率歷史示意圖Fig.1 Diagram of reactor power history

對式(6)在時間0～TN內積分可得：

(7)

分別在圖1所示的N個時間分段內考慮式(7)右端積分，可得：

exp[-λj(TN-T1)]+

exp[-λj(TN-Tn)]+NN+

(8)

其中：M為所考慮的裂變產物或錒系核素種類，系統程序處理中一般可將裂變產物或錒系核素根據其衰變常量的大小歸并為若干組，此時M即表示組數；N1，N2，N3，…，NN分別為T1，T2，T3，…，TN時刻的裂變功率。根據式(8)可求得N1，N2，N3，…，NN，最終對于N段功率歷史的情形，第j種核素的衰變熱的初始值為：

λjHj(TN)=EjN1(1-e-λjT1)e-λj(TN-T1)+

e-λj(TN-Tn)+EjNN(1-e-λj(TN-TN-1))

(9)

總的衰變熱Pd為各組產物衰變熱之和：

(10)

最后，初始化過程中如果某一功率歷史分段末尾已有的衰變熱高于該時段的給定功率，則取該時段末尾裂變功率為0。例如，如果反應堆的功率歷史中有一段時間處于零功率停堆狀態，則初始化計算中該時間段末尾對應的裂變功率取值為0。

2 模型驗證

模型驗證分為兩部分，第一部分選取ANSI/ANS-5.1—2005標準[14]中的數據，不考慮裂變功率影響，只計算停堆后的衰變熱，驗證不同功率歷史下衰變熱的計算；第二部分選取中國實驗快堆(CEFR)的設計數據，考慮裂變功率的影響，驗證衰變熱的計算。

2.1 基于ANSI/ANS-5.1—2005標準的驗證

美國核學會(ANS)衰變熱計算標準(ANSI/ANS-5.1—2005)是ANS于2005年發布的針對輕水堆的衰變熱計算標準，該標準給出了235U、238U、239Pu和241Pu等核素發生單次裂變及輻照足夠長時間(1013s)后衰變熱隨時間的變化[14]，標準中給出的數據可用于驗證所開發的模型。

1) 單一功率歷史驗證

考慮某一核素輻照足夠長時間(1013s)后衰變熱的變化，可近似等效為以該種核素為堆芯的反應堆以固定功率運行1013s后停堆，計算停堆后衰變熱隨時間的變化，功率歷史如圖2所示。

只有一段功率歷史時，t時刻的衰變熱為：

(11)

其中，N1為Top時刻的裂變功率，由下式計算：

(12)

計算時裂變產物分組參考標準劃分為23組。每組產物對應的衰變常量λj、無量綱數Ej可由標準給出的數據求得。FR-Sdaso程序計算結果和ANS標準數據的比較列于表1、2，因標準中只給出了4位有效數字，所以程序計算結果也取4位有效數字。

表1 輻照1013 s后停堆FR-Sdaso計算和ANS標準的相對衰變熱Table 1 Relative decay heat calculated by FR-Sdaso and in ANS standard after 1013 s irradiation shutdown

綜上可見，所開發衰變熱模型的計算結果與ANSI/ANS-5.1—2005標準結果相比，在所計算的時長范圍內(0～109s)，采用科學計數法保留4位有效數字時，最大相對偏差約為0.1%，結果符合良好，初步證明了所開發模型的正確性。

2) 多段功率歷史驗證

SAS4A/SASYS-1程序[3]是美國阿貢國家實驗室開發的鈉冷快堆系統分析程序，該程序已經過較為充分的驗證，廣泛應用于鈉冷快堆系統分析計算[15-16]。為進一步驗證所開發模型對功率歷史影響模擬的正確性，以ANSI/ANS-5.1—2005標準中衰變熱計算采用的衰變常量λj和無量綱數Ej為輸入條件，假設了一種反應堆功率歷史，分別計算了235U、238U、239Pu和241Pu經歷所假設的功率歷史后衰變熱的變化，并與SAS4A/SASYS-1的計算結果[17]進行比較。

假設功率歷史如圖3所示，劃分為4段，反應堆首先滿功率運行13周，停堆維持零功率1周，然后維持50%功率1周，最后維持滿功率1 d后停堆，計算此后衰變熱隨時間的變化。4種核素的計算結果及偏差列于表3、4。可看出，FR-Sdaso和SAS4A/SASYS-1計算的最大相對偏差約為5.6×10-8，結果符合良好，初步證明了所開發模型的正確性。

表2 輻照1013 s后停堆FR-Sdaso與ANS標準相對衰變熱的相對偏差Table 2 Relative deviation between FR-Sdaso and ANS after 1013 s irradiation shutdown

圖3 多段功率歷史示意圖Fig.3 Schematic of multi period power history

表3 多段功率歷史下FR-Sdaso和SAS4A/SASYS-1計算所得相對衰變熱Table 3 Relative decay heat calculated by FR-Sdaso and SAS4A/SASYS-1 under multi period power history

表4 多段功率歷史下FR-Sdaso與SAS4A/SASYS-1計算結果相對偏差Table 4 Relative deviation between FR-Sdaso and SAS4A/SASYS-1 results under multi period power history

2.2 CEFR停堆算例驗證

CEFR停堆算例考慮裂變功率對衰變熱的影響，采用CEFR的設計參數，分別采用FR-Sdaso和SAS4A/SASYS-1程序計算了CEFR緊急停堆過程中衰變熱的變化，功率歷史考慮反應堆滿功率連續運行80 d后緊急停堆，主要的輸入參數列于表5～7。

表5 CEFR停堆算例主要輸入參數Table 5 Main input parameter of CEFR shutdown calculation case

表6 緩發中子的衰變常量和有效緩發中子份額Table 6 Decay constant and effective fraction of delayed neutron

為單獨驗證功率計算，不考慮反饋反應性，緊急停堆過程中控制棒引入的反應性采用式(13)模擬。

表7 各裂變產物或錒系核素的衰變常量和功率份額Table 7 Decay constant and power fraction of each group of fission product or actinide

(13)

其中：ρ(t)為緊急停堆后t時刻引入的反應性，pcm；t為緊急停堆后的時刻，s。

圖4為FR-Sdaso和SAS4A/SASYS-1計算所得裂變功率和衰變熱。可看出，FR-Sdaso和SAS4A/SASYS-1程序計算所得裂變功率的最大相對偏差在10-6量級，衰變熱的最大相對偏差在10-8量級，結果符合良好。衰變熱的最大相對偏差出現在緊急停堆初期，控制棒下插到底后，相對偏差出現波動，后續隨著裂變功率和衰變熱的衰減，相對偏差逐漸減小，計算結

圖4 FR-Sdaso和SAS4A/SASYS-1計算所得裂變功率和衰變熱Fig.4 Fission power and decay heat calculated by FR-Sdaso and SAS4A/SASYS-1

果初步證明了所開發模型的正確性。

3 裂變功率對衰變熱的影響分析

為進一步分析裂變功率在緊急停堆瞬態過程中對衰變熱的影響，采用FR-Sdaso計算了不考慮裂變功率影響時的衰變熱，對比結果如圖5所示。圖中衰變熱Y和衰變熱N分別表示考慮和不考慮裂變功率影響的結果。緊急停堆過程中，裂變功率的衰減需要一定的時間，在此期間仍不斷有裂變產物或錒系核素產生，因此考慮裂變功率對衰變熱的影響后，計算所得的衰變熱更高，兩者的最大相對偏差(相對偏差=(衰變熱N-衰變熱Y)/衰變熱Y)出現在停堆初期，約為-7%，出現最大偏差的時刻被低估的衰變熱約為0.29 MW，后續隨緊急停堆后裂變功率和衰變熱的下降，相對偏差逐漸減小，約100 s后相對偏差下降至-1%以下。

圖5 停堆初期裂變功率對衰變熱的影響Fig.5 Influence of fission power on decay heat during initial shutdown

4 結論

本文開發了一種可用于鈉冷快堆系統分析程序的衰變熱計算模型，該模型可考慮裂變功率和功率歷史對衰變熱的影響。并通過與ANSI/ANS-5.1—2005標準和SAS4A/SASYS-1程序計算結果對比，初步驗證了所開發模型的正確性。基于中國實驗快堆的設計數據，分析了緊急停堆過程中裂變功率的變化對衰變熱的影響，結果表明計算停堆初期衰變熱時應考慮裂變功率的影響。

本模型的開發為自主化池式鈉冷快堆系統分析程序FR-Sdaso提供了衰變熱模型，后續將繼續對模型進行分析和驗證，為我國大型鈉冷快堆的設計分析提供技術支持。