鉛基快堆氧化物燃料性能分析程序FUTURE開發與初步驗證

王 冠，顧 龍,3，于 銳，王 挺，王 兆，袁 和，惲 迪,*

(1.中國科學院 近代物理研究所，甘肅 蘭州 730000；2.中國科學院大學，北京 100049；3.蘭州大學，甘肅 蘭州 730000；4.西安交通大學，陜西 西安 710049)

21世紀以來，第4代核能系統國際論壇(GIF)面向未來能源需求，從安全性、經濟性、可持續性等多個角度出發，選擇了6種最具前景的先進反應堆系統[1]。目前，以鉛冷快堆、鈉冷快堆等為代表的第4代核能技術正逐步走向工程化應用的軌道，這種趨勢對快堆的理論認識、計算模擬和安全設計提出了很高的要求。相比于金屬燃料，U-Pu混合氧化物(MOX)燃料的穩定性較好，在高燃耗時腫脹更小，并具有較高的熔點，因此MOX燃料是當前快中子增殖堆的首選燃料[2-3]。

在反應堆內部高溫、強輻照環境下，燃料內部發生著從微觀原子尺度到宏觀結構尺度的演化，同時也經歷著微秒量級的中子學反饋和數年的服役壽期[4]。相比于熱中子堆，快堆內中子能譜更硬、功率密度和燃耗深度更大，造成了MOX燃料元件內部存在顯著的燃料重構、元素遷移、裂變氣體釋放、芯包相互作用以及各種化學效應等[5-6]。這些現象之間強烈的非線性特征和復雜的耦合關系使得燃料元件性能分析只能通過計算機數值模擬實現。因此從20世紀開始，許多國家相繼開發了一系列針對液態金屬快堆氧化物燃料元件性能分析程序從而實現對燃料元件性能和行為的數值模擬和定量分析，例如歐洲超鈾元素研究所(TUI)開發的TRANSURANUS[7]，美國麻省理工學院(MIT)開發的FEAST-OXIDE[8]，法國原子能與可替代能源委員會(CEA)開發的GERMINAL[9]以及日本原子能委員會開發的CEDAR[10]等。

為實現對鉛基快堆以及CiADS[11-12]次臨界堆氧化物燃料元件性能和行為的模擬，中國科學院近代物理研究所反應堆室和西安交通大學NUSOL團隊聯合開發了鉛基快堆氧化物燃料棒性能分析程序FUTURE。該程序基于串行的半隱式耦合求解方法，可實現對鉛基快堆氧化物燃料的熱工-力學-輻照-化學等多物理場耦合分析計算。相比于國際上的1.5維快堆燃料性能分析程序，一方面FUTURE程序采用了兩步分析法耦合多物理場計算，在實現較高精度求解的同時提升了計算效率，具有后發優勢；另一方面，國際上的燃料性能分析程序多針對鈉冷快堆進行開發和研究，能夠實現鉛基快堆燃料性能分析的程序較少，而FUTURE程序初期即面向鉛基快堆進行開發，并添加了目前較為先進的液態鉛鉍(LBE)腐蝕模型。

本文針對FUTURE程序中實現的快譜條件下氧化物燃料特異性行為進行模型的初步驗證，為后續基于程序實現鉛基快堆氧化物燃料棒性能分析和模擬計算奠定基礎。

1 程序介紹

1.1 計算建模

燃料元件是反應堆內具有獨立結構的燃料使用單元，也是燃料組件的核心構件。典型的液態金屬快堆燃料元件呈棒狀結構，外部為細長的金屬包殼管(長約為2～3 m，直徑約為5～10 mm，厚度約為0.4～0.6 mm[6])，內部包容有燃料活性區、氣腔及其他具有定位、測量、屏蔽等功能的部件。FUTURE程序將燃料元件的燃料區和氣腔作為主要模擬計算區域，采用1.5維軸向堆疊假設和有限體積法進行簡化和區域離散，如圖1所示。離散后的控制體為相互嵌套的節點環，從而生成計算網格。為了對快堆MOX燃料棒的復雜行為進行模擬計算，程序基于Fortran95語言進行模塊化編程，各計算模塊實現相應的計算和分析功能并進行數據傳遞，如圖2所示。

圖1 FUTURE程序燃料計算模型Fig.1 Fuel calculation model of FUTURE code

圖2 FUTURE程序模塊框架Fig.2 Framework of FUTURE code

1.2 計算流程

在快堆高溫、大溫度梯度和強中子輻照作用下，MOX燃料內部存在顯著的孔隙遷移、燃料重構、氧和钚的徑向重分布、輻照腫脹、蠕變、裂變氣體釋放和芯包相互作用等現象。這些現象以及它們之間的相互作用十分復雜，為了對燃料在服役期間的各種行為進行模擬和定量計算，需要通過燃料元件性能分析程序對燃料、包殼和冷卻劑的重要特性隨輻照時間的演變加以追蹤，從而預測燃料棒在各種運行工況下的性能[13-14]。基于串行的半隱式計算架構，程序采用燃料棒全域熱力分析和局部行為模擬耦合的兩步分析法，能夠較好地反映和處理燃料棒在堆內服役過程中的多物理場耦合和多尺度演化。FUTURE程序計算流程如圖3所示。

圖3 FUTURE程序計算流程Fig.3 Flow chart of FUTURE code

2 數學物理模型

2.1 傳熱分析

燃料棒服役期間的許多行為都強烈地依賴燃料棒內部的溫度分布并且能夠反作用于傳熱過程，為精確計算燃料棒內部的溫度分布和演化過程，程序通過耦合力學和燃料行為模型，采用高斯-賽德爾迭代方法對燃料棒傳熱方程進行數值求解。

芯塊和包殼內部的熱量傳遞過程由一維固體熱傳導方程描述：

(1)

其中：ρ為材料密度，kg/m3；c為材料比熱容，J/(kg·K)；T為溫度，K；t為時間，s；r為徑向位置，m；λ為熱導率，W/(m·K)；Q為燃料釋熱率，W/m3。

選擇圍繞單根燃料棒的獨立冷卻劑通道作為流體研究區域，將冷卻劑流道類比成內壁面加熱、外壁面絕熱的環管，采用單通道分析計算冷卻劑軸向傳熱過程，冷卻劑軸向能量方程為：

(2)

其中：Al為燃料棒冷卻劑流通面積，m2；w為冷卻劑軸向流速，m/s；z為軸向位置，m；qf為包殼-冷卻劑界面熱通量，W/m。

2.2 力學分析

在力學模擬計算中，主要考慮了材料在外部載荷及高溫輻照環境下發生的彈性、熱膨脹、腫脹、重定位、密實化、蠕變和塑性變形等。為了降低程序力學計算的復雜度，計算引入如下基本假設[15]：1) 連續性假設；2) 軸對稱假設；3) 廣義平面應變假設；4) 各向同性假設；5) 小變形假設；6) 準靜態假設。

力學計算從基本方程出發，采用位移法建立應力、應變與位移的對應關系，并通過全主元高斯-若當消去法進行數值求解，力學控制方程如下。

1) 徑向平衡方程

該方程由極坐標下微元體的受力平衡得到，不考慮徑向體積力時方程為：

(3)

其中，σr、σθ為徑向和環向主應力，Pa。

2) 幾何方程

基于小變形假設和廣義平面應變假設，柱坐標系下燃料棒模型適用的幾何方程為：

(4)

其中：εr、εθ、εz為徑向、環向和軸向的主應變；u、w為徑向、軸向位移，m；C為常數。

3) 物理本構方程

采用廣義胡克定律描述應力-應變的本構關系，在彈性項的基礎上添加了熱膨脹、腫脹、密實化、重定位以及塑性和蠕變應變，其形式為：

(5)

2.3 燃料重構

在高溫和大溫度梯度作用下，快堆內運行的氧化物燃料由于初始孔隙的徑向遷移和原始晶粒的不斷長大，其幾何形態和微觀結構都會發生很大變化。程序通過求解孔隙率守恒方程[5]計算和模擬燃料內部孔隙的遷移和中心開孔現象，并耦合傳熱和力學分析，復現快堆柱狀芯塊向環形芯塊的演化過程。孔隙率守恒方程為：

(6)

其中：P為燃料孔隙率；vp為孔隙遷移速率，m/s。

程序采用Sens[16]的孔隙遷移速率公式：

vp=5×10-16×(0.988+6.395×10-6T+

3.543×10-9T2+3×10-12T3)·

(7)

其中：ΔH為蒸發焓，ergs/mol；p0為蒸汽壓，dynes/cm2；R為理想氣體常數。

2.4 元素重分布

MOX燃料中鈾和钚存在許多價態，所以這些元素偏離準確的化學計量是允許的，這種偏離會對燃料芯塊的性能產生很大影響。程序基于固體熱擴散的遷移動力學理論，計算了MOX燃料中氧和钚的重分布過程，其中用于間隙氧原子和氧空位的遷移動力學方程[17]為：

(8)

其中：cO為間隙氧原子/氧空位的原子分數；DO為擴散系數，m2/s；QO為有效摩爾輸運熱，J/mol。

用于钚熱擴散的遷移動力學方程[18]為：

(9)

其中：cPu為MOX燃料分子中钚元素的原子分數；DU-Pu為鈾钚互擴散系數，m2/s；QU-Pu為有效摩爾輸運熱，J/mol。

2.5 LBE腐蝕

為評估LBE對燃料元件包殼的影響，FUTURE程序基于Zhang[19-20]開發的長周期氧化-腐蝕動力學模型，計算和預測不銹鋼包殼在服役期間的腐蝕和氧化情況。

在外部Fe3O4磁鐵層未完全去除之前，尖晶石層并不與LBE接觸，所以此時尖晶石層的生長符合拋物線定律：

(10)

其中：δSp為尖晶石層厚度，m；kp為氧化常數，m2/s。

外部Fe3O4磁鐵層的生長動力學方程為：

(11)

其中：δFe3O4為磁鐵層厚度，m；ρFe3O4、ρSt、ρSp、ρLBE為Fe3O4磁鐵層、不銹鋼基體、尖晶石層、LBE冷卻劑密度，kg/m3；FFe,Fe3O4、FFe，St、FFe，Sp為Fe3O4磁鐵層、不銹鋼基體、尖晶石層中Fe的質量分數；Rm為Fe的質量傳輸率，m/s。

當LBE流動造成的Fe去除速率小于通過尖晶石層Fe的擴散速率時，LBE對尖晶石層無腐蝕作用，此時尖晶石層遵循線性生長律：

(12)

其中，t0為外層Fe3O4磁鐵層去除時間，s。

隨著尖晶石層按照線性律不斷生長，其厚度最終會趨于一個臨界厚度，在此厚度下LBE流動造成的Fe去除速率等于通過尖晶石層Fe的擴散速率，此時LBE對尖晶石層開始存在腐蝕作用，因此尖晶石的生長開始受到Tedmon模型的控制：

(13)

其中，RSp為尖晶石層腐蝕速率，m/s。

3 初步驗證

鉛鉍快堆氧化物燃料元件性能分析程序FUTURE包含有傳熱、力學、燃料重構、元素重分布、LBE腐蝕等主要計算模塊。由于各模塊內部計算過程復雜，模塊間耦合接口眾多，為了測試各計算模塊的穩定性和計算精度，保證程序全耦合條件下的穩定運行和正確計算，分別通過與文獻算例、基準公式、現有程序進行對比，開展了程序各主要計算模塊的初步驗證。

3.1 傳熱驗證

傳熱模塊的驗證采用與FEAST-OXIDE溫度驗證相同的基準題[21]，該基準題針對UO2+He2+HT9燃料元件傳熱材料體系，基于簡單的工況和物性假設，模擬計算特定幾何結構的燃料元件內部溫度分布。該基準題的溫度分布滿足下列解析表達式：

(14)

其中：T(r)為溫度分布，K；q?為體功率密度，W/m3；ro為燃料芯塊外徑，m；k為燃料熱導率，W/(m·K)；Ts為燃料芯塊外表面溫度，K。

FUTURE程序傳熱分析對比驗證如圖4所示。圖4結果表明，在快堆定常工況條件下，燃料內部溫度呈現拋物線分布，當線功率密度較高時，燃料中心溫度逼近熔點。基于溫度基準題線功率66 kW/m和70 kW/m兩個工況的計算結果分析，程序的傳熱計算比較準確，能較好地模擬燃料元件內部的傳熱過程和溫度分布。

圖4 FUTURE程序傳熱分析對比驗證Fig.4 Comparison of FUTURE code heat analysis

3.2 力學驗證

3.2.1芯塊非均勻熱彈性驗證 通過FUTURE程序和商用的多物理場耦合平臺COMSOL對預設溫度場的柱狀芯塊和環形芯塊非均勻熱彈性算例進行對比，驗證程序基于廣義平面應變假設和1.5維軸向堆疊模型力學核心算法的可靠度和準確性。力學驗證算例的輸入參數列于表1，材料物性熱膨脹系數α、楊氏模量E、泊松比ν和密度ρ如下式所列。

α=1.0×10-5K-1

(15)

E=2.334×1011×(1-1.091 5×

10-4T)×(1.0-2.752×

(1-0.925))×(1+0.15×0.2)

(16)

ν=(1-0.2)×(0.316+(T-300.0)/

2 800.0×(0.5-0.316))+0.2×

(0.276+(T-300.0)/2 800.0×(0.5-0.276))

(17)

(18)

表1 力學驗證算例輸入參數Table 1 Input parameter of mechanical analysis case

柱狀芯塊和環形芯塊非均勻熱彈性算例預設溫度分布如圖5所示，力學計算對比如圖6所示。由圖6可見：在高溫和大溫度梯度條件下，FUTURE程序力學計算與COMSOL的計算結果符合較好，其中反映燃料變形的位移計算基本無偏差；von Mises應力計算趨勢相同，結果相近，FUTURE程序柱狀芯塊和環形芯塊應力計算與COMSOL計算結果的最大相對偏差分別為8.27%和6.48%。基于上述燃料非均勻熱彈性算例的計算結果分析，FUTURE程序針對不同結構的靜力學計算具有較高的準確度和可靠性，能實現快堆MOX燃料幾何結構從柱狀向環形演化過程的力學分析和模擬。

a——柱狀芯塊；b——環形芯塊

a——柱狀芯塊；b——環形芯塊

3.2.2芯包機械相互作用驗證 在快堆高溫、高燃耗條件下，燃料元件內部存在著比較嚴重的芯包機械相互作用(PCMI)。通過FUTURE程序和COMSOL對環形芯塊芯包接觸工況算例進行對照，驗證程序執行PCMI計算的準確度和可靠性。環形芯塊PCMI工況算例溫度分布如圖7所示，力學計算對比如圖8所示。

圖7 環形芯塊PCMI工況溫度分布Fig.7 Given temperature distribution of annular pellet PCMI case

由圖8可見，在接觸工況條件下環形燃料的最大應力出現在芯塊的內壁面和外壁面處，而接觸壓作為芯塊和包殼外部載荷影響了內部的應力分布和結構變形。對于芯塊而言，環向應力和軸向應力在內壁面呈壓應力，外壁面處呈拉應力狀態。這種應力狀態會導致芯塊外緣出程序計算的芯包內部各應力分量和位移與COMSOL結果符合較好，其中徑向應力計算基本無偏差，環向應力和軸向應力最大相對偏差為8.19%和11.11%。基于接觸工況的對比驗證，表明FUTURE程序具備處理和模擬芯包機械相互作用的能力。

圖8 環形芯塊PCMI工況力學計算對比Fig.8 Comparison of annular pellet PCMI analysis

圖9、10為FUTURE程序加載典型快堆功率譜并開啟熱力耦合(熱彈性)功能計算得到的UO2燃料中心線溫度和接觸壓與COMSOL模擬的結果對比。由圖9可見，FUTURE程序計算得到的燃料中心線溫度隨功率的演化與COMSOL計算結果趨勢完全相同，最大相對偏差為3.9%。由圖10可見，FUTURE程序計算的接觸壓與COMSOL的相比變化趨勢相同，數值相近，但仍存在一定偏差，最大相對偏差在20%左右。但FUTURE程序計算的接觸后芯包邊界的貫穿深度比COMSOL計算的更好，其貫穿深度更小。從方法和應用上來看，COMSOL是基于有限元法開發的三維多物理場耦合平臺，在力學計算上具有更高的精度和靈活性，而FUTURE程序采用了簡化的有限差分方法處理力學問題，受方法和網格量限制，精度無法比及COMSOL，但在現有框架下已能較好地模擬快堆燃料棒的熱力耦合問題。

圖9 快堆工況下UO2燃料中心線溫度演化對比Fig. 9 Comparison of UO2 fuel central temperature for fast reactor

圖10 快堆工況下UO2燃料FCMI演化對比Fig.10 Comparison of UO2 fuel FCMI simulation for fast reactor

3.3 燃料重構驗證

FUTURE程序燃料重構計算選用Sens[16]開發的氧化物燃料孔隙遷移速率公式，并與文獻[16]模擬結果進行對比驗證，結果如圖11所示。

由圖11可見，快堆氧化物燃料內部初始分布均勻的孔隙會發生沿徑向的定向遷移，這種遷移過程不僅使孔隙率重新分布，影響燃料局部的熱力學特性，也會使燃料中心在短時間內富集較多的孔隙導致中心孔的出現。圖11a、b結果表明，程序計算的孔隙率遷移與文獻結果符合較好，能較好地反映孔隙率的演化過程。圖11c顯示在孔隙率演化逐漸達到穩態的情況下，FUTURE程序計算的由于孔隙率遷移導致的燃料靠近中心處密實化更大，分析其原因可能為：1) 由于文獻[16]中計算模擬的輸入條件不明確，程序的孔隙遷移計算采用弱對照的方法，擬合并匹配輸入條件，可能會導致計算出現一定的偏差；2) 程序孔隙遷移采用有限體積法進行計算，守恒性好，而文獻[16]中采用的數值方法和離散格式未知，可能導致計算結果精度有所差異。

a——30 min；b——1 h；c——2 h

3.4 元素重分布驗證

3.4.1氧的重分布驗證 FUTURE程序中氧的重分布計算選用Sari等[17]開發的氧重布遷移動力學模型和計算關系式，并與文獻[17]模擬結果進行對比，用以驗證程序氧分布計算的準確度和可靠性。

快堆氧化物燃料中發生的氧的重分布現象影響了O/M比(氧元素與重金屬元素的含量之比)的變化，而O/M比的變化將會直接改變燃料的各種特性并引發一系列現象。快堆氧化物燃料O/M比演化計算對比如圖12所示。圖12中，R為外半徑。由圖12可見：在快堆工況下，對于超化學計量的UO2或MOX燃料，氧的定向遷移會使燃料熱端的含氧量增加，同時冷端含氧量相應減少；在欠化學計量燃料中，氧的重布則以反向發生，使冷端濃集氧，熱端失去氧。通過FUTURE程序計算的不同初始O/M比條件下氧的重分布結果與文獻[17]結果對比發現，程序中超化學計量氧的重分布計算與文獻[17]符合較好，但在欠化學計量燃料(O/M<2.0)的計算中，推測文獻中可能采用固定O/M=2.0為外邊界條件，而FUTURE程序采用第二類邊界條件，導致外側區域O/M比的計算有所偏離。

紅色線為FUTURE程序計算結果，黑色線為文獻[17]結果

3.4.2钚的重分布驗證 為模擬和分析快堆MOX燃料中钚元素的重分布現象，FUTURE程序中钚的重分布計算選用Bober等[18]開發的基于固體熱擴散機理的動力學模型和經驗公式，并與TRANSURANUS程序的模擬結果[22]進行對比，結果如圖13所示。

圖13 快堆MOX燃料钚元素遷移計算對比Fig.13 Comparison of MOX fuel Pu migration calculation for fast reactor

在快堆MOX燃料運行過程中，由于高溫和大溫度梯度的作用會在燃料內部發生鈾、钚成分的離析。圖13表明，在快堆啟堆的幾十天內，燃料中心處钚的濃度大幅攀升，這種钚的重分布過程會造成局部燃料成分的變化，不僅影響燃料熱物理性能，也會由于钚在中心處的富集導致燃料中心溫度大幅上升。基于Bober等[18]的熱擴散動力學模型，FUTURE程序模擬計算結果與TRANSURANUS程序的符合較好，表明能較好處理和模擬快堆MOX燃料钚的重分布過程。

3.5 LBE腐蝕驗證

圖14為 FUTURE程序計算的高流速(v=5 m/s)條件下長周期LBE腐蝕情況與文獻[20]模擬結果的對比。由圖14可見，程序腐蝕計算判定準確，能夠較好地切換不同計算模型，從而模擬和追蹤長周期HT-9和T-91不銹鋼包殼雙層氧化層的生長與腐蝕情況。由于文獻[20]中部分輸入條件不明確，計算采用的數值方法未知，導致程序模擬計算結果在較長周期的腐蝕演化計算中不斷出現累積偏差，10 a腐蝕總量的相對偏差為8.82%(HT-9)和9.09%(T-91)。考慮到快堆燃料元件服役壽期遠小于10 a，且程序計算具有一定的保守性，因此認為程序LBE腐蝕計算能對不銹鋼包殼的腐蝕情況進行較好地模擬和評估。

圖14 LBE腐蝕計算對比Fig.14 Comparison of LBE corrosion calculation

4 結論

本文基于串行的半隱式耦合求解方法開發了鉛基快堆氧化物燃料性能分析程序FUTURE。FUTURE程序采用兩步分析法，能夠實現鉛基快堆氧化物燃料棒全域熱力分析與局部行為模型的多物理場耦合計算。通過與文獻算例、基準公式和現有程序的對照，初步驗證了程序對快堆氧化物燃料各分離效應的模擬方法，以及對各計算模塊求解的準確度、穩定性和誤差來源進行了分析。結果表明，FUTURE程序能模擬快堆穩態工況條件下燃料元件內部的溫度演化、結構變形、應力分布和相互作用，并實現燃耗模擬、燃料重構、組分遷移、裂變氣體釋放和LBE腐蝕等計算內容。其中熱力分析、燃料重構、氧和钚元素的遷移等計算具有較高的計算精度，燃料服役期內的LBE腐蝕計算也具有較好的計算結果。除此以外，FUTURE程序在實現快堆燃料熱力耦合計算中能夠保證一定的計算精度和穩定性，為接下來針對程序的全耦合驗證和工程應用提供可靠的模擬計算內容，同時也為快堆穩態工況條件下氧化物燃料的行為演化和性能分析奠定基礎。