快堆系統分析程序FASYS堆芯分析模塊驗證

王 晉，張東輝

(中國原子能科學研究院，北京 102413)

快堆系統分析程序可分析快堆電廠在各種瞬態過程中系統的響應[1]，進行快堆研究的國家和機構開展了大量的快堆系統分析程序開發和驗證。目前快堆系統分析程序主要有兩大類：一類是由水堆系統分析程序經過模型升級和修改成為可適用于快堆分析的系統分析程序，如RELAP5-3D?[2]、TRACE[3]、CATHARE[4]、ATHLET[5]等；另一類是專門為快堆系統分析而研發的系統分析程序，如SSC-L[6]、SSC-P[7]、SAS4A/SASSYS-1[8]、SAM[9]、DINROS[10]、OASIS[11]、SSC-K[12]、MARS-LMR[13]等。

FASYS程序是中國原子能科學研究院自主開發的快堆系統分析程序，包含堆芯分析模塊、一二回路模塊、事故余熱排出系統模塊等，已對FASYS程序的中子物理模型、水力模型、熱工模型和整體模型進行了初步驗證[14-15]。目前FASYS程序正用于中國示范快堆反應性的意外變化類事故的分析，主要關注堆芯關鍵現象的分析計算。

DINROS程序(俄羅斯快堆系統分析程序)曾用于БН-600和БН-800快中子反應堆的安全分析報告的編寫及БН-350反應堆安全保護手段的研制，該程序已在БН-350、БН-600、БP-10等裝置或反應堆上經過驗證計算。中國實驗快堆曾引進DINROS程序進行安全分析報告的編制。SAS4A/SASSYS-1程序(以下簡稱SAS程序)由美國阿貢國家實驗室開發，其中包括SAS4A程序和SASSYS-1程序。SAS4A程序是嚴重事故程序，SASSYS-1程序是快堆系統分析程序，可用于分析與系統相關的典型事故的序列，也是EBR-Ⅱ仿真機的計算引擎，被廣泛地用于快堆事故的分析中[8]。

FASYS程序堆芯分析模塊包括點堆、衰變熱、反應性反饋、堆芯通道熱工水力模型等，本文通過點堆方程解析解算例、DINROS程序超功率算例、SAS4A/SASSYS-1程序點堆與衰變熱算例、反應性反饋算例、燃料棒與冷卻劑換熱算例，進行FASYS程序堆芯分析模塊的所有關鍵模型的驗證，并對堆芯分析模塊的計算偏差進行初步評估。

1 FASYS程序堆芯分析模塊介紹

FASYS程序堆芯分析模塊包括點堆、衰變熱、反應性反饋、堆芯通道熱工水力模型等。FASYS程序采用3階Hermite插值多項式法[16]求解點堆動態方程。在計算反應堆功率時，除考慮裂變功率外，還考慮裂變產物和錒系元素的衰變功率，并考慮了反應堆連續運行史對上述衰變功率的影響。FASYS程序堆芯通道熱工水力模型采用單通道模型，堆芯可劃分為任意數目的通道來模擬燃料組件、屏蔽組件、不銹鋼組件等。文獻[14]給出了FASYS程序點堆、衰變熱、堆芯通道熱工水力模型的具體方程。FASYS程序可計算的反應性分項包括輸入卡定義的引入反應性、控制棒引入的反應性、燃料多普勒反應性反饋、燃料和包殼軸向膨脹反應性反饋、冷卻劑密度變化反應性反饋、堆芯徑向膨脹反應性反饋。其中，燃料多普勒反應性反饋、燃料和包殼軸向膨脹反應性反饋、冷卻劑密度變化反應性反饋提供集總參數模型和空間分布模型兩個選項。

2 驗證結果

2.1 點堆方程解析解算例

本算例將采用點堆方程解析解對FASYS程序的點堆模型進行驗證，解析解的來源為公開發表的文獻[16]，包括3種不同情況下解析解與FASYS程序數值解的比較。

情況1：

M=1，β=0.006 5，
λ=0.08 s-1，Λ=1.0×10-8s

式中：M為緩發中子組數；β為緩發中子有效份額；λ為先驅核衰變常量；Λ為瞬發中子壽命；ρ為引入的外部反應性。

情況2：

M=15，β=0.007 330 19，
Λ=2.15×10-4s(動態參數見文獻[16])

情況3：

M=1，β=0.006 5，
λ=0.08 s-1，Λ=1.0×10-8s

ρ(t)=0.003sin(πt/10)

情況1和情況2均為0 s時刻引入0.002階躍反應性，但中子動態參數不同，情況1僅有1組緩發中子，情況2有15組緩發中子。表1列出情況1解析解與FASYS程序計算值(N(t)/N(0))的對比，共計算時間步長分別為0.1、0.01、0.001、0.000 1和0.000 01 s等5種情況，由于0 s時刻引入的階躍反應性較大，采用的時間步長不能滿足程序設定的精度要求，程序中自動減小了時間步長，以上5種情況在0 s時刻的實際計算時間步長均為10-6s量級，0 s時刻之后采用原設置的時間步長。從表1可知，當時間步長為0.1 s時，相對偏差為10-3量級，當時間步長減小1個量級，計算精度約提高1個量級。表2列出情況2解析解與FASYS程序計算值的對比，可看出當時間步長為0.001 s時，相對偏差為10-5量級，時間步長越小，計算精度越高。情況3引入隨時間變化的正弦反應性，中子動態參數與情況1相同。表3列出情況3解析解與FASYS程序計算值的對比，可看出當時間步長為0.01 s時，FASYS程序計算值即可與解析解保持6位有效數字的一致。

通過情況1、2、3的計算對比可發現，FASYS程序采用的3階Hermite插值法計算的數值結果在一般情況下較為精確，對于引入較大階躍反應性的情況需采取較小的時間步長才能保證計算精度，在實際計算中應評估引入的反應性并進行時間步長影響分析，選取合適的時間步長以保證計算精度。

表1 不同時間步長下情況1解析解與FASYS程序計算值的對比Table 1 Comparison between analytical solution and FASYS code calculated value of case 1 at different time steps

表2 不同時間步長下情況2解析解與FASYS程序計算值的對比Table 2 Comparison between analytical solution and FASYS code calculated value of case 2 at different time steps

表3 情況3解析解與FASYS程序計算值的對比Table 3 Comparison between analytical solution and FASYS code calculated value of case 3

2.2 DINROS程序超功率算例

本算例將采用中國實驗快堆[10]“在堆各種狀態下補償棒非規定位移”預計運行事件分析中對0.1%額定功率運行工況的計算結果進行對比驗證，原分析結果由DINROS程序計算。工況假設如下：1) 相對流量取0.25；2) 堆芯冷卻劑入口溫度為250 ℃；3) 堆內無反饋；4) 棒價值取平衡態初期的值，引入反應性7.77×10-4Δk/k，在5 s內線性引入。

表4列出了FASYS程序計算的主要事故序列與原事故序列的對比，可看出，FASYS程序計算的達到堆相對功率保護參數和最大值的時間與DINROS程序計算值一致，堆相對功率最大值的相對偏差為8.9×10-4。

表4 主要事故序列對比Table 4 Comparison of main accident sequences

圖1 反應堆相對功率的結果對比Fig.1 Comparison of normalized power results

圖1、2分別為反應堆相對功率和反應堆周期的FASYS程序計算值與DINROS程序計算值的對比，從結果可看出，FASYS程序計算值與DINROS程序計算值符合很好。

2.3 SAS程序點堆與衰變熱計算算例

本算例將采用SAS程序作為校驗程序對FASYS程序的點堆模型和衰變熱模型進行驗證，計算建?；谥袊鴮嶒灴於?，假設反應堆已滿功率連續運行80 d，0～6 s向反應堆線性引入負反應性，反應性速率為-728.68 pcm/s，6 s時負反應性引入結束，共引入負反應性-4 372.08 pcm。不考慮反應性反饋，計算停堆過程中的裂變功率和衰變功率。

圖2 反應堆周期的結果對比Fig.2 Comparison of nuclear reactor period results

表5列出了衰變功率相對值、總功率相對值，可看出，衰變功率的最大相對偏差出現在0 s時刻，為7.4×10-7。總功率的最大相對偏差出現在前期，主要來源于兩個程序求解點堆方程的方法不同，后期由于衰變功率占總功率的絕大部分，總功率的相對偏差也逐步減小?？傮w而言，本算例中FASYS程序與SAS程序的衰變功率最大相對偏差為10-7量級，而總功率最大相對偏差為10-6量級。

表5 衰變功率相對值、總功率相對值Table 5 Normalized decay power results and normalized total power results

2.4 SAS程序反應性反饋算例

本算例將采用SAS程序作為校驗程序對FASYS程序的反應性反饋模型進行驗證，將驗證具有空間分布模型的燃料多普勒反應性反饋、燃料和包殼軸向膨脹反應性反饋、冷卻劑密度變化反應性反饋。本算例中FASYS程序與SAS程序采用空間分布相同的鈉空泡反應性系數、多普勒常數、燃料質量增加引入的反應性、包殼質量增加引入的反應性。

情況1堆芯零功率，其初始溫度均為358 ℃，0～50 s時間內堆芯入口溫度線性上升至458 ℃，50～100 s堆芯入口溫度一直保持在458 ℃；情況2堆芯零功率，堆芯初始溫度均為358 ℃，0～50 s時間內堆芯入口溫度線性下降至258 ℃，50～100 s堆芯入口溫度一直保持在258 ℃。

表6列出了情況1在100 s時反應性反饋值的FASYS程序計算值與SAS程序計算值的對比。其中，鈉密度變化引入的反應性反饋計算值相對偏差絕對值為1.3×10-3，軸向膨脹引入的反應性反饋計算值相對偏差絕對值為6.6×10-4，多普勒效應引入的反應性反饋計算值相對偏差絕對值為1.1×10-3。表7列出了算例2在100 s時反應性反饋值的FASYS程序計算值與SAS程序計算值的對比。其中，鈉密度變化引入的反應性反饋計算值相對偏差絕對值為1.3×10-3，軸向膨脹引入的反應性反饋計算值相對偏差絕對值為1.1×10-3，多普勒效應引入的反應性反饋計算值相對偏差絕對值為1.1×10-3。總體來說，FASYS程序計算的反應性反饋與SAS程序相比，其相對偏差為10-3量級。

表6 情況1反應性反饋計算值對比Table 6 Comparison of reactivity feedback results for case 1

表7 情況2反應性反饋計算值對比Table 7 Comparison of reactivity feedback results for case 2

2.5 SAS程序燃料棒與冷卻劑換熱計算算例

本算例將采用SAS程序作為校驗程序對FASYS程序的堆芯通道熱工水力模型進行驗證，采用相同的燃料組件幾何參數、節點劃分、功率、流量、軸向功率分布、包殼物性、燃料物性。其中，燃料棒軸向共劃分24個節點，芯塊徑向等距劃分為5個節點，包殼徑向劃分為3個節點，冷卻劑徑向為1個節點。對100%功率穩態情況下和100%功率線性升功率到110%功率瞬態情況下SAS程序與FASYS程序計算的燃料溫度、包殼溫度、冷卻劑溫度進行結果對比。

圖3為100%功率穩態情況下SAS程序與FASYS程序計算的燃料最高溫度、包殼中壁溫度、冷卻劑溫度的軸向分布結果對比。其中，燃料最高溫度、包殼中壁溫度、冷卻劑溫度的最大偏差絕對值均為0.2 ℃。

圖3 組件軸向溫度的穩態結果對比Fig.3 Comparison of subassembly axial temperature steady state results

瞬態計算中，0 s時刻燃料組件為100%功率然后線性升功率，8 s時刻達到110%功率后保持不變，期間流量保持不變，時間步長為0.001 s。圖4為FASYS程序與SAS程序的瞬態燃料最高溫度結果對比，瞬態情況下最大偏差絕對值為0.1 ℃。圖5為FASYS程序與SAS程序的瞬態包殼中壁最高溫度與冷卻劑最高溫度結果對比，瞬態情況下包殼中壁最高溫度最大偏差絕對值為0.4 ℃，冷卻劑最高溫度最大偏差絕對值為0.7 ℃。

圖4 瞬態燃料最高溫度的結果對比Fig.4 Comparison of fuel maximum temperature transient results

圖5 瞬態包殼中壁最高溫度和冷卻劑最高溫度的結果對比Fig.5 Comparison of cladding middle wall maximum temperature and coolant maximum temperature transient results

分析時間步長對圖4、5瞬態溫度計算的影響，分別計算時間步長為0.01、0.001和0.000 1 s等3種情況，如表8所列，時間步長0.01 s的計算結果與時間步長0.000 1 s相比溫度最大偏差為0.07 ℃，時間步長0.001 s的計算結果與時間步長0.000 1 s相比溫度最大偏差為0.01 ℃。

表8 不同時間步長下瞬態溫度偏差對比Table 8 Comparison of transient temperature deviations at different time steps

注：溫度偏差指與時間步長為0.000 1 s的計算結果的偏差

分析網格劃分對燃料棒溫度計算的影響，由于軸向節點劃分方式來自物理專業給出的軸向功率分布，此處僅分析芯塊徑向節點劃分數目對計算的影響。分別計算芯塊徑向節點數目為5、6、9、11、13、15、17、21、31、41等10種情況，對比前9種情況的燃料最高溫度計算值與第10種情況計算值的偏差，如圖6所示，可看出隨著徑向節點數目的增加，燃料最高溫度計算值減小，偏差亦減小，當徑向節點大于21時，燃料溫度計算偏差小于0.1 ℃，基本實現徑向網格無關性。

圖6 不同徑向節點數的燃料溫度偏差對比Fig.6 Deviation comparison of fuel temperature results with different numbers of radial nodes

分析功率和流量的偏差對圖4和圖5瞬態溫度計算的影響，分別計算瞬態最大功率增加1%、增加0.1%、瞬態流量減小1%、減小0.1%等情況，如表9所列，功率增加1%和流量減小1%引起的包殼溫度、冷卻劑溫度偏差相當，而功率增加1%引起的燃料溫度偏差為13.5 ℃，遠大于流量減小1%引起的燃料溫度偏差1.0 ℃。

表9 功率、流量偏差對瞬態溫度的影響Table 9 Effect of power and flow deviations on transient temperature

注：溫度偏差指與SAS程序原瞬態計算結果的偏差

3 堆芯分析模塊的計算偏差評估

以FASYS程序進行中國實驗快堆的調節棒非規定位移事故分析為例，對堆芯分析模塊的整體計算偏差進行分析。首先給出調節棒非規定位移事故分析的事故描述和事故假設，反應堆處于正常運行期間，1根調節棒從底部失控提升到頂，1根調節棒價值為0.001 748 Δk/k，假設反應性線性引入，反應性引入速率為11.66 pcm/s，反應性引入時間從0 s開始到15 s結束。事故假設初始狀態100%額定功率上疊加2.5%，堆芯冷卻劑入口溫度在額定值360 ℃上疊加3 ℃，緊急停堆共引入負反應性-0.016 8 Δk/k，停堆反應性在1.4 s內引入完成。表10列出計算時間步長取0.001 s時的主要事故序列。

表10 主要事故序列Table 10 Major accident sequences

分析計算時間步長h對結果的影響，圖7a為時間步長取0.01、0.001 s的功率計算結果與時間步長0.000 1 s的功率計算結果的相對偏差，可看出在7.0 s時，由于緊急停堆短時間內(1.4 s)引入較大的負反應性，計算偏差急劇上升，時間步長取0.01 s的計算偏差從1×10-4瞬間上升到4×10-2，時間步長取0.001 s的計算偏差從1×10-5瞬間上升到1×10-3，緊急停堆結束后，由于不再引入較大的負反應性，時間步長取0.01 s的計算偏差降低到約1×10-3，時間步長取0.001 s的計算偏差降低到約1×10-4。圖7b為時間步長取0.01、0.001 s的燃料最高溫度計算結果與時間步長0.000 1 s的燃料最高溫度計算結果的偏差，時間步長取0.01 s的燃料最高溫度計算值偏差約10 ℃，時間步長取0.001 s的燃料最高溫度計算值偏差約0.2 ℃，由表8、9和圖6可知，不同時間步長的燃料最高溫度計算值偏差主要來自功率的計算偏差。圖7c為時間步長取0.01、0.001 s的包殼、冷卻劑最高溫度計算結果與時間步長0.000 1 s的包殼、冷卻劑最高溫度計算結果的偏差，時間步長取0.01 s的包殼最高溫度計算值偏差約1 ℃、冷卻劑最高溫度計算值偏差約0.4 ℃，時間步長取0.001 s的包殼最高溫度計算值偏差約0.1 ℃、冷卻劑最高溫度計算值偏差約0.1 ℃，不同時間步長的包殼、冷卻劑最高溫度計算值偏差應來自燃料溫度的計算偏差。

綜合第2章對堆芯分析模塊各模型的驗證結果及本章結果可知，在引入較大反應性的情況下時間步長對點堆方程的求解精度有較大影響，對溫度計算偏差的影響最大，徑向節點的劃分數目對燃料溫度計算的影響次之，隨著徑向節點數目的增加，燃料最高溫度計算值減小。在此類事故的分析中，應首先評估引入的外部反應性，然后進行時間獨立性分析，選取合適的時間步長，進行網格無關性分析，得到滿足計算偏差要求的徑向節點數目。

圖7 不同時間步長下功率、燃料最高溫度以及包殼、冷卻劑最高溫度偏差對比Fig.7 Comparison of deviations for power, fuel maximum temperature, and cladding and coolant temperatures at different time steps

4 結論

通過點堆方程解析解算例、DINROS程序超功率算例、SAS4A/SASSYS-1程序點堆與衰變熱算例、反應性反饋算例、燃料棒與冷卻劑換熱算例，對FASYS程序的堆芯分析模塊的所有關鍵模型進行驗證，結果均符合良好。對堆芯分析模塊各關鍵模型和整體計算偏差進行了初步評估，在引入較大反應性的情況下點堆方程的求解精度依賴于時間步長，功率計算值的偏差直接影響溫度計算的準確性，徑向節點的劃分數目對燃料溫度計算的影響次之，隨著徑向節點數目的增加，燃料最高溫度計算值減小，為進行中國示范快堆反應性的意外變化類事故分析和計算偏差估計提供參考。