基于MOOSE平臺(tái)的棒狀燃料元件性能分析程序開發(fā)與驗(yàn)證

鄧超群，向烽瑞，賀亞男，牛鈺航，巫英偉，田文喜，秋穗正，蘇光輝

(西安交通大學(xué) 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

福島核電站事故后，事故容錯(cuò)燃料(ATF)成為研究熱點(diǎn)[1]，涂層Zr包殼[2]、復(fù)合SiC包殼[3]和U3Si2芯塊[4]等多種ATF設(shè)計(jì)被提出。傳統(tǒng)燃料性能分析程序FRAPCON[5]、FRAPTRAN[6]和TRANSURANUS[7]等難以實(shí)現(xiàn)ATF復(fù)雜結(jié)構(gòu)(如多層包殼結(jié)構(gòu))的建模和分析。因此，亟需開發(fā)新的燃料性能分析程序?qū)TF進(jìn)行性能分析及設(shè)計(jì)改進(jìn)。

當(dāng)前，美國(guó)愛達(dá)荷國(guó)家實(shí)驗(yàn)室(INL)開發(fā)了多物理場(chǎng)燃料性能分析程序BISON[8]，并針對(duì)UO2-Zr燃料開展了大量驗(yàn)證，但由于技術(shù)保密等原因，我國(guó)難以獲取其使用權(quán)限。Deng和He等[9-10]分別基于FRAPCON4.0、FRAPTRAN2.0程序添加了多層SiC穩(wěn)態(tài)、瞬態(tài)分析模型，實(shí)現(xiàn)了多層SiC包殼燃料棒全堆芯穩(wěn)態(tài)及典型棒瞬態(tài)分析，但程序采用差分離散算法，計(jì)算結(jié)果并不精細(xì)。

為實(shí)現(xiàn)對(duì)ATF行為的模擬，西安交通大學(xué)核反應(yīng)堆熱工水力研究室(NuTHeL)基于開源的多物理場(chǎng)有限元平臺(tái)MOOSE[11]開發(fā)棒狀燃料元件性能分析程序BEEs[4]，可實(shí)現(xiàn)二維柱坐標(biāo)下的穩(wěn)態(tài)棒狀燃料性能分析。由于國(guó)際上ATF實(shí)驗(yàn)數(shù)據(jù)較少，目前主要針對(duì)BEEs程序中添加的UO2-Zr燃料物性與行為模型開展驗(yàn)證工作，為后續(xù)基于BEEs程序?qū)崿F(xiàn)ATF性能分析奠定基礎(chǔ)。

1 控制方程

BEEs程序中燃料熱傳導(dǎo)方程如下：

(1)

變形計(jì)算中假定芯塊、包殼均處于靜態(tài)平衡狀態(tài)：

(2)

其中：σ為柯西應(yīng)力張量，Pa；f為體積力，N·kg-1。

2 數(shù)學(xué)物理模型

當(dāng)前BEEs程序中添加的UO2-Zr燃料物性及行為模型主要來自MATPRO報(bào)告[12]、FRAPCON[5]和BISON理論手冊(cè)[13]。

2.1 燃耗計(jì)算

燃耗計(jì)算采用TUBRNP模型[14]，局部平均核子密度為：

(3)

2.2 間隙換熱

燃料包殼和芯塊的間隙換熱模型[15]如下：

hgap=hg+hs+hr

(4)

其中：hgap為間隙換熱系數(shù)，W·m-2·K-1；hg為氣體換熱系數(shù)，W·m-2·K-1；hs為芯塊包殼接觸換熱系數(shù)，W·m-2·K-1；hr為芯塊外表面、包殼內(nèi)表面輻射換熱系數(shù)，W·m-2·K-1。

2.3 包殼氧化腐蝕

燃料包殼氧化層厚度計(jì)算采用以下模型[16]：

(5)

其中：i和i+1分別為上一時(shí)間步和當(dāng)前時(shí)間步；s為包殼氧化層厚度，m；Δw為氧化物增重，g·cm-2；t為時(shí)間，d；A為氧化系數(shù)，A=6.3×109μm3·d；R為氣體常數(shù)；stran為模型轉(zhuǎn)換厚度，stran=2.0×10-6m；Ti為氧化層交界面溫度，Ti=To+q″s/λ，K；To為冷卻劑溫度，K；q″為包殼外表面熱通量，W·m-2；λ為ZrO2導(dǎo)熱系數(shù)，W·m-1·K-1；γ為氧化層增重與厚度的轉(zhuǎn)換因子，γ=0.678 9 cm-3·g；k0=11 863+3.5×104×(1.91×10-15φ)0.24g·cm-2·d，φ為快中子通量密度；Q1和Q2為不同階段包殼氧化激活能，cal·mol-1。對(duì)于M5包殼，Q1=27 446 cal·mol-1，Q2=29 816 cal·mol-1；對(duì)于ZIRLO包殼，Q1=27 446 cal·mol-1，Q2=27 354 cal·mol-1；對(duì)于Zr-2和Zr-4包殼，Q1=32 289 cal·mol-1，Q2=27 354 cal·mol-1。

2.4 密實(shí)化、重定位與輻照腫脹

UO2燃料芯塊密實(shí)化、重定位和輻照腫脹采用ESCORE[17]和MATPRO[12]經(jīng)驗(yàn)?zāi)Ｐ停?/p>

(6)

(7)

exp(-0.016 2(2 800-T))·

exp(-0.017 8ρBu)

(8)

(9)

(10)

(11)

其中：T為燃料溫度，℃；q′為平均線功率，kW·ft-1。

2.5 芯塊和包殼蠕變

UO2芯塊蠕變采用MATPRO FCREEP[12]模型，包殼熱蠕變和輻照蠕變計(jì)算分別采用Power-law模型[18]和Hoppe模型[19]：

(12)

(13)

(14)

2.6 裂變氣體釋放

BEEs采用Forsberg-Massih模型[20]模擬裂變氣體釋放過程，裂變氣體在晶粒內(nèi)的擴(kuò)散方程如下：

(15)

其中：C為晶粒內(nèi)裂變氣體濃度，m-3；Deff為裂變氣體有效擴(kuò)散系數(shù)，m2·s-1；β為裂變氣體產(chǎn)生率，m-3·s-1。

邊界條件為：

(16)

其中：r為距晶粒中心的距離，m；a為晶粒半徑，m；b為重溶率，%；λ為重溶半徑，m；Nb為晶界氣體濃度，m-2。

當(dāng)晶界氣體濃度達(dá)到飽和值時(shí)就會(huì)釋放。晶界氣體飽和濃度Ns為：

(17)

其中：γ為表面張力，J·m-2；θ為二面角半角，(°)；F為氣泡體積轉(zhuǎn)換因子；Vc為晶界臨界覆蓋率，%；KB為玻爾茲曼常量；pext為氣泡外壓，Pa。

3 程序框架

MOOSE是基于有限元方法、JFNK算法的開源平臺(tái)。當(dāng)前MOOSE平臺(tái)中已開發(fā)了多種物理場(chǎng)和豐富的擴(kuò)展接口，支持多場(chǎng)耦合和全耦合、全隱性求解計(jì)算。

BEEs程序總體框架如圖1所示。通過繼承MOOSE平臺(tái)C++類模板，本文基于MOOSE平臺(tái)的張量力學(xué)模塊、熱傳導(dǎo)模塊和接觸模塊，添加了壓水堆燃料相關(guān)行為模型和非線性燃料元件熱膨脹系數(shù)、導(dǎo)熱系數(shù)等物性模型，實(shí)現(xiàn)了棒狀燃料元件性能分析程序的開發(fā)。求解計(jì)算時(shí)，采用有限元庫(kù)Libmesh完成網(wǎng)格讀取與有限元方程離散，由PETSc進(jìn)行串行或并行求解。

圖1 BEEs程序框架Fig.1 Diagram of BEEs code

4 程序模塊驗(yàn)證與校核

在前期開發(fā)工作中，通過對(duì)比BISON和BAQUS程序，驗(yàn)證了BEEs程序針對(duì)UO2-Zr燃料溫度、空腔內(nèi)壓、間隙寬度變化的計(jì)算功能[4]。本文通過對(duì)比FRAPCON程序和實(shí)驗(yàn)數(shù)據(jù)，進(jìn)一步驗(yàn)證BEEs針對(duì)長(zhǎng)期穩(wěn)態(tài)工況下燃料燃耗、氧化層厚度和燃料溫度模擬的準(zhǔn)確性。

4.1 燃耗計(jì)算校核

隨燃耗的加深，徑向功率分布形狀由于邊緣效應(yīng)將發(fā)生顯著變化，而燃料性能參數(shù)和裂變產(chǎn)物腫脹等行為模擬都與局部燃耗密切相關(guān)。FRAPCON程序[5]的燃耗計(jì)算模塊經(jīng)大量驗(yàn)證，計(jì)算結(jié)果可靠。本文針對(duì)同一基準(zhǔn)題進(jìn)行模擬，計(jì)算UO2燃料徑向功率分布，計(jì)算參數(shù)列于表1。圖2示出徑向功率分布因子對(duì)比。由圖2可見，在燃耗初始階段及較高燃耗階段，BEEs程序計(jì)算結(jié)果與FRAPCON的均符合較好，校核了程序燃耗模塊的正確性。

表1 燃耗基準(zhǔn)題計(jì)算參數(shù)Table 1 Calculation parameter of burnup benchmark problem

4.2 氧化層計(jì)算驗(yàn)證

Oconee Rod 15309實(shí)驗(yàn)[21]在Oconee PWR中完成了5個(gè)輻照周期，并對(duì)燃耗末期的氧化層厚度進(jìn)行了測(cè)量，實(shí)驗(yàn)中燃料棒設(shè)計(jì)參數(shù)列于表2，燃料類型為UO2，包殼類型為Zr-4。圖3示出氧化層厚度軸向分布對(duì)比。由圖3可知，BEEs程序計(jì)算結(jié)果與FRAPCON程序計(jì)算結(jié)果和實(shí)驗(yàn)值均符合良好，驗(yàn)證了BEEs程序計(jì)算氧化層厚度結(jié)果的合理性。

a——零燃耗；b——平均燃耗為33 MW·d/kgU圖2 徑向功率分布因子對(duì)比Fig.2 Comparison of radial power factor

表2 Oconee Rod 15309實(shí)驗(yàn)設(shè)計(jì)參數(shù)Table 2 Design parameter of Oconee Rod 15309 test

圖3 氧化層厚度軸向分布對(duì)比Fig.3 Comparison of axial distribution of oxide thickness

4.3 燃料溫度計(jì)算驗(yàn)證

Halden IFA432實(shí)驗(yàn)[22]是于Halden反應(yīng)堆進(jìn)行的針對(duì)BWR-6燃料棒的長(zhǎng)期穩(wěn)態(tài)輻照實(shí)驗(yàn)，實(shí)驗(yàn)中由下端伸入的熱電偶測(cè)得燃料棒中心溫度，實(shí)驗(yàn)中燃料棒Rod1和Rod3的設(shè)計(jì)參數(shù)列于表3，燃料類型為UO2，包殼類型為Zr-2。圖4示出Rod1和Rod3的燃料中心溫度的變化。由圖4可知，BEEs程序計(jì)算結(jié)果與實(shí)驗(yàn)值總體符合良好，與FRAPCON計(jì)算結(jié)果基本一致，驗(yàn)證了BEEs程序在長(zhǎng)期穩(wěn)態(tài)工況下燃料溫度變化預(yù)測(cè)的準(zhǔn)確性。

表3 Halden IFA432實(shí)驗(yàn)燃料棒Rod1和Rod3的設(shè)計(jì)參數(shù)Table 3 Design parameter of Rod1 and Rod3 for Halden IFA432 test

圖4 Rod1(a)和Rod3(b)燃料棒中心溫度對(duì)比Fig.4 Fuel rod centerline temperature comparison of Rod1 (a) and Rod3 (b)

5 BR3 Rod實(shí)驗(yàn)算例整體分析

BR3 Rod實(shí)驗(yàn)[23]是于比利時(shí)BR3反應(yīng)堆開展的高燃耗輻照實(shí)驗(yàn)。本文采用BR3 Rod實(shí)驗(yàn)進(jìn)行整體分析并與FRAPCON及實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，進(jìn)一步驗(yàn)證BEEs程序?qū)核寻魻钊剂闲袨槟M的準(zhǔn)確性。實(shí)驗(yàn)中燃料棒的相關(guān)參數(shù)列于表4，燃料類型為UO2，包殼類型為Zr-4，平均線功率的變化如圖5所示。

表4 BR3 Rod實(shí)驗(yàn)燃料棒設(shè)計(jì)參數(shù)Table 4 Design parameter of fuel rod for BR3 Rod test

圖5 燃料平均線功率隨時(shí)間的變化Fig.5 Average linear power of fuel vs. time

模擬時(shí)將碟形加倒角燃料芯塊層疊形狀簡(jiǎn)化單一燃料棒，燃料棒頂部為平端面。網(wǎng)格劃分時(shí)在燃料徑向劃分20個(gè)節(jié)點(diǎn)，軸向劃分100個(gè)節(jié)點(diǎn)，包殼徑向劃分4個(gè)節(jié)點(diǎn)，軸向劃分80個(gè)節(jié)點(diǎn)，網(wǎng)格類型為QUAD4。

圖6示出包殼內(nèi)表面、芯塊外表面和芯塊中心最高溫度隨時(shí)間的變化。結(jié)合圖5燃料平均線功率變化可知，在279 d時(shí)，由于功率激增，燃料棒整體迅速升溫，在熱膨脹作用下芯塊、包殼間隙減小，直至最終間隙閉合，如圖7所示。

圖6 包殼內(nèi)表面、芯塊外表面和芯塊中心最高溫度隨時(shí)間的變化Fig.6 Maximum temperature at cladding inner surface, fuel outer surface and fuel centerline vs. time

圖7 間隙尺寸隨時(shí)間變化Fig.7 Gap width vs. time

選取燃料棒中心高度處的包殼內(nèi)表面和芯塊中心線溫度進(jìn)行對(duì)比，結(jié)果分別如圖8所示。由圖8可知，BEEs程序預(yù)測(cè)結(jié)果合理，其計(jì)算結(jié)果與FRAPCON的基本一致。

裂變氣體釋放計(jì)算結(jié)果對(duì)比示于圖9。輻照前期，產(chǎn)生的少量裂變氣體擴(kuò)散至晶界沉積。在279 d時(shí)，由于功率激增，裂變氣體生成量顯著增加，同時(shí)燃料棒溫度提升到較高水平，使得裂變氣體擴(kuò)散率較大，積累在晶界的氣體濃度超過飽和值，裂變氣體開始釋放。BEEs程序預(yù)測(cè)裂變氣體釋放份額與FRAPCON程序的符合良好，但總體高于實(shí)驗(yàn)測(cè)量值，計(jì)算結(jié)果偏保守。

圖8 包殼內(nèi)表面(a)和芯塊(b)在中心高度處溫度的對(duì)比Fig.8 Comparison of temperature of cladding inner surface (a) and fuel centerline (b) at axial mid-plane

圖9 裂變氣體釋放份額對(duì)比Fig.9 Comparison of fission gas release portion

6 總結(jié)

本文基于開源的多物理場(chǎng)有限元平臺(tái)MOOSE開發(fā)了棒狀燃料元件性能分析程序BEEs。通過與FRAPCON程序和相關(guān)實(shí)驗(yàn)測(cè)量值的對(duì)比，初步驗(yàn)證了程序的穩(wěn)態(tài)工況模擬能力。結(jié)果表明，BEEs能實(shí)現(xiàn)燃料燃耗、氧化層厚度和燃料溫度的準(zhǔn)確計(jì)算及裂變氣體釋放的合理預(yù)測(cè)，且能針對(duì)典型壓水堆燃料元件在長(zhǎng)期穩(wěn)態(tài)工況下的整體性能進(jìn)行合理分析。