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

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

(西安交通大學(xué) 動力工程多相流國家重點實驗室，陜西 西安 710049)

福島核電站事故后，提出了眾多事故容錯燃料(ATF)包殼和芯塊選型[1]，旨在提高燃料抵抗事故的能力。國內(nèi)外學(xué)者針對瞬態(tài)事故下ATF行為開展了廣泛研究。部分研究對FRAPTRAN、TRANSURANUS等傳統(tǒng)燃料性能分析程序進行改造，分析了U3Si2[2]、FeCrAl[3-4]和SiC[5]等在瞬態(tài)事故工況下的性能，但難以實現(xiàn)如陶瓷基包覆顆粒彌散燃料(FCM)[6]、涂層包殼[7]等復(fù)雜結(jié)構(gòu)燃料的精細建模。因此需開發(fā)新的瞬態(tài)燃料性能分析程序，準(zhǔn)確評估瞬態(tài)事故條件下ATF的性能。

目前，基于商用有限元平臺如COMSOL[6]、ABAQUS[8]、ADINA[9]等開發(fā)的復(fù)雜結(jié)構(gòu)ATF性能分析程序均缺乏瞬態(tài)事故分析功能。美國愛達荷國家實驗室(INL)基于更具擴展性的開源多物理場有限元平臺MOOSE[10]開發(fā)的BISON程序[11]能夠模擬瞬態(tài)事故下ATF燃料行為，但國內(nèi)無法獲得其使用權(quán)。

西安交通大學(xué)核反應(yīng)堆熱工水力研究室(NuTHeL)基于MOOSE平臺開發(fā)了棒狀燃料元件穩(wěn)態(tài)性能分析程序BEEs[12]，該程序包含了燃耗、冷卻劑焓升、芯塊密實化、重定位、熱膨脹、晶粒長大、裂變氣體釋放、裂變產(chǎn)物腫脹、間隙換熱、燃料內(nèi)壓演化等行為模型及非線性燃料元件熱膨脹系數(shù)、導(dǎo)熱系數(shù)等物性模型。

為實現(xiàn)反應(yīng)性引入事故(RIA)和冷卻劑喪失事故(LOCA)瞬態(tài)事故工況下的燃料性能精確模擬，本文針對UO2-Zr燃料棒，基于棒狀燃料元件穩(wěn)態(tài)性能分析程序BEEs進行瞬態(tài)功能擴展，并開展程序模塊驗證與整體功能驗證。

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

BEEs程序基于MOOSE平臺的張量力學(xué)模塊、接觸模塊和熱傳導(dǎo)模塊進行開發(fā)，底層的網(wǎng)格離散和求解計算依賴于PETSc、Hypre、MPICH和LibMesh等功能庫，其整體框架示于圖1。為實現(xiàn)RIA和LOCA瞬態(tài)事故下燃料性能分析，在BEEs程序中添加了包殼彈塑性、高溫蠕變、高溫相變、高溫氧化和包殼失效等模型。

圖1 BEEs程序結(jié)構(gòu)框架Fig.1 Diagram of BEEs code framework

1.1 彈塑性模型[13]

在發(fā)生屈服前，用Hooke定律描述包殼變形：

σ=ε·E

(1)

其中：σ為應(yīng)力，MPa；E為彈性模量，MPa；ε為應(yīng)變。

當(dāng)包殼的應(yīng)力達到屈服極限后，用冪次定律描述包殼的塑性變形：

(2)

1.2 高溫蠕變模型

發(fā)生LOCA時，包殼在燃料棒升溫和內(nèi)壓作用下發(fā)生高溫蠕變，導(dǎo)致局部產(chǎn)生較大變形。BEEs程序中針對包殼的高溫蠕變速率計算以Norton方程[14]的形式給出：

(3)

表1 Zr4包殼蠕變常數(shù)Table 1 Creep constant of Zr4 cladding

1.3 高溫相變模型

隨著包殼溫度上升，鋯合金晶體從六角α相結(jié)構(gòu)轉(zhuǎn)變?yōu)榱⒎叫桅孪嘟Y(jié)構(gòu)。鋯包殼的相變速率可用下式[15]表示：

(4)

其中：y為β相鋯的體積分數(shù)，%，0≤y≤1；t為時間，s；k為速率常數(shù)，s-1；ys為平衡態(tài)β相鋯的體積分數(shù)，%。模型適用于包殼升溫速率小于100 K·s-1。

(5)

(6)

其中：w為局部氫濃度，ppm；k0為動力常數(shù)，s-1；E為等效激活能，J；kb為玻爾茲曼常數(shù)，J·K-1；km為常數(shù)，s-1。對于Zr4包殼，k0=60 457+18 129|Z|(Z為加熱(或降溫)速率，K·s-1，0.1≤|Z|≤100 K·s-1)；E/kb=16 650 K。

常數(shù)km取值為：

(7)

相變轉(zhuǎn)換起始溫度為：

(8)

(9)

1.4 高溫氧化模型

LOCA下燃料元件溫度快速上升，包殼氧化速率遠大于正常穩(wěn)態(tài)工況，其氧化動力學(xué)方程以拋物線形式給出[16]：

(10)

(11)

其中：δ為氧化物厚度，m；w為氧化物質(zhì)量，kg；Aδ和Aw為氧化系數(shù)，m2·s-1；Qδ和Qw為激活能，J·mol-1；R為通用氣體常數(shù)，J·mol-1·K-1；Ti為金屬-氧化層界面溫度，K，Ti=To+q″δ/λ，To為冷卻劑溫度，K，q″為包殼外表面熱通量，W·m-2，λ為氧化鋯的導(dǎo)熱系數(shù)，W·m-1·K-1。

高溫氧化模型適用于Ti>673 K。模型常數(shù)取值列于表2。當(dāng)673 K1 900 K時，采用Prater-Courtright關(guān)系式；當(dāng)1 800 K

表2 氧化關(guān)系式參數(shù)Table 2 Parameter of oxidation relationship

1.5 包殼失效準(zhǔn)則[14]

包殼在RIA中由于芯塊和包殼產(chǎn)生的強烈相互作用(PCMI)產(chǎn)生的較大應(yīng)力導(dǎo)致失效，在LOCA中由于發(fā)生塑性變形使其鼓脹破裂。BEEs程序中采用包殼局部應(yīng)力限制準(zhǔn)則和塑性“失穩(wěn)”準(zhǔn)則來判定包殼失效。事故過程中，包殼狀態(tài)滿足其一即判定包殼失效。

1) 包殼局部應(yīng)力限制準(zhǔn)則

(12)

其中：σθ為包殼環(huán)向應(yīng)力，MPa；σb為包殼破裂臨界應(yīng)力，MPa；a和b為模型常數(shù)，根據(jù)表3中參數(shù)插值獲得；Δη為包殼氧化物質(zhì)量分數(shù)增量，%；rcl，o和rcl，i分別為包殼外徑和內(nèi)徑，m；w為局部氫濃度，ppm；ρZr為包殼密度，kg·m-3；rmet，o為未氧化包殼厚度，m。

表3 局部應(yīng)力限制準(zhǔn)則參數(shù)Table 3 Parameter of local stress limiting criterion

2) 塑性“失穩(wěn)”準(zhǔn)則

(13)

2 瞬態(tài)模型驗證

為確保程序中模型添加的正確性，本文首先對BEEs程序添加的彈塑性計算、高溫蠕變、高溫相變和高溫氧化模塊進行單一模塊驗證。

2.1 包殼彈塑性求解驗證

彈塑性模塊驗證中，將邊長為1 m的Zr樣品保持恒溫，下方固定，上方施加恒定速率的位移，計算參數(shù)列于表4，獲得拉伸過程中的總應(yīng)變與應(yīng)力關(guān)系，并與實驗數(shù)據(jù)[13]進行比較，如圖2所示。在達到屈服極限后樣品硬化，應(yīng)力應(yīng)變曲線斜率發(fā)生變化，達到抗拉極限后樣品斷裂，程序計算終止。由圖2可知，BEEs程序預(yù)測值與文獻計算值和實驗結(jié)果符合較好，由于理論模型未對Zr合金屈服后的均勻塑性變形和強化階段進行區(qū)分，因此計算結(jié)果和實驗曲線的趨勢有所不同，在均勻塑性變形和強化階段相比于實驗值的最大相對誤差分別為4.7%和11.3%。

表4 應(yīng)力應(yīng)變關(guān)系式驗證算例參數(shù)Table 4 Parameter of validation case of stress-strain correlation

2.2 高溫蠕變模型驗證

Hardy tube實驗[17]是為評估Zr4包殼在假定LOCA實驗下的鼓脹行為而進行的瞬態(tài)實驗。實驗測量了升溫速率為25 K·s-1時不同內(nèi)壓下包殼的變形情況。圖3示出1.4 MPa內(nèi)壓條件下環(huán)向應(yīng)變的對比，計算參數(shù)列于表5，包殼軸向和徑向分別劃分20和4個節(jié)點。當(dāng)包殼處于安全限值內(nèi)(環(huán)向應(yīng)變低于2.78%)時，BEEs程序在環(huán)向應(yīng)變低于1%范圍內(nèi)計算相對誤差較大，最大相對誤差為26.5%，當(dāng)環(huán)向應(yīng)變高于1%時最大相對誤差為8.6%。當(dāng)超過包殼安全限值時相對誤差增大，但此時已判定包殼失效，不影響燃料性能的有效評估。需要注意的是，由于BEEs程序在700～900 K之間進行了低溫?zé)崛渥僛18-19]與高溫蠕變的插值，而BISON程序均采用高溫蠕變計算關(guān)系式，因此在該溫度區(qū)間內(nèi)計算結(jié)果比BISON程序略高。

圖3 內(nèi)壓1.4 MPa時包殼環(huán)向應(yīng)變對比Fig.3 Comparison of cladding strain with internal pressure at 1.4 MPa

表5 Hardy tube實驗設(shè)計參數(shù)Table 5 Design parameter of Hardy tube test

2.3 包殼相變模型驗證

實驗[15]中將Zr合金樣品置于膨脹儀中升溫(或降溫)，同時進行樣品膨脹(收縮)量的數(shù)據(jù)采集(精度±0.1 μm)獲得熱膨脹量曲線，再通過杠桿法或切線法計算β相轉(zhuǎn)化比例。表6列出包殼相變驗證算例參數(shù)。圖4示出升溫(降溫)速率為10 K·s-1時，BEEs程序計算的加熱時β相體積分數(shù)及平衡態(tài)β相體積分數(shù)隨溫度的變化，并與實驗值進行對比，計算結(jié)果合理。圖5對比了0.1～100 K·s-1升溫速率范圍內(nèi)β相體積分數(shù)到達10%的溫度(T10%β)，BEEs程序計算值相對誤差均在1.85%以內(nèi)。

表6 包殼相變驗證算例參數(shù)Table 6 Parameter of validation case of cladding phase transformation

圖4 β相體積分數(shù)隨溫度的變化Fig.4 Comparison of volume fraction of β phase as a function of temperature

圖5 不同升溫速率下T10%β的對比Fig.5 Comparison of T10%β as a function of heating rate

2.4 包殼高溫氧化模型驗證

Yoo等[20]針對Zr的氧化動力學(xué)規(guī)律進行了實驗研究。將Zr合金樣品經(jīng)過打磨拋光和酸溶液去銹后，置于石墨恒溫爐(精度為±1 K)中，用熱天平(精度為±1 μg)測量氧化物的增重情況，并將氧化物增重轉(zhuǎn)換成氧化層厚度的增長，實驗中氣體溫度工況范圍為773～1 073 K。圖6示出BEEs程序計算的氧化層厚度與實驗值[20-21]的對比。由圖6可知，高溫氧化模型適用于在溫度較高(大于973 K)、氧化層較厚(大于10 μm)的情況，最大相對誤差為11.4%。當(dāng)溫度較低且氧化層厚度較小時，計算誤差增大。針對低溫范圍內(nèi)的穩(wěn)態(tài)工況計算，BEEs程序采用同F(xiàn)RAPCON一致的LEISTIKOW模型[22]。在溫度為973 K時，計算值與實驗值偏離可能是由于采用氧化模型未考慮摻雜元素的影響或該組實驗過程中出現(xiàn)氧化層脫落的情況。

圖6 氧化層厚度變化對比Fig.6 Comparison of oxide thickness

3 整體驗證分析

3.1 RIA模擬驗證

CABRI-RIA Na-2實驗[23]是IRSN同EDF、Framatome、CEA和NRC在法國合作開展的RIA瞬態(tài)實驗。本文基于該算例對BEEs程序的RIA整體分析功能進行評估。實驗中燃料棒的相關(guān)參數(shù)列于表7，燃料類型為UO2，包殼類型為Zr4。實驗功率脈沖歷史如圖7所示，RIA功率脈沖峰值出現(xiàn)在0.08 s，平均功率峰值約為25 000 kW·m-1。模擬時在芯塊軸向劃分340個節(jié)點，包殼徑向劃分4個節(jié)點，軸向劃分350個節(jié)點，網(wǎng)格類型為QUAD4，如圖8所示。

圖7 實驗功率脈沖歷史Fig.7 History of experimental power pulse

圖8 Na-2實驗算例網(wǎng)格劃分Fig.8 Mesh partition of case Na-2

表7 Na-2實驗設(shè)計參數(shù)Table 7 Design parameter of Na-2 test

Na-2實驗基礎(chǔ)輻照工況末期燃料平均燃耗約為25 MW·d/kgU，為準(zhǔn)確模擬高燃耗下芯塊顯著的邊緣效應(yīng)，在其徑向劃分不同網(wǎng)格節(jié)點，對比功率峰值節(jié)點(PPN)處芯塊溫度徑向分布，如圖9所示。芯塊徑向劃分35個節(jié)點時與24節(jié)點和40節(jié)點的計算最大相對誤差分別為0.4%和0.1%，滿足計算精度要求。

圖9 不同節(jié)點劃分情況下PPN處芯塊溫度徑向分布對比Fig.9 Radial distribution of pellet temperature at PPN under different nodes partitions

圖10示出PPN處芯塊中心溫度和包殼內(nèi)表面溫度的變化。其中BISON程序和FRAPTRAN程序計算結(jié)果來自文獻[5]。模擬時BEEs程序采用與FRAPTRAN程序相同的冷卻劑邊界，因此兩者溫度變化趨勢總體一致。由圖10可見，在RIA初期，BEEs和FRAPTRAN程序的溫度計算值符合良好，后期溫度的差異可能與程序間包殼變形程度計算不同有關(guān)。

圖10 PPN處芯塊中心溫度及包殼內(nèi)表面溫度預(yù)測對比Fig.10 Comparison of pellet centerline and cladding inner surface temperatures at PPN

在RIA過程中PCMI使包殼發(fā)生較大塑性變形，長期冷卻后仍有殘余變形。圖11示出冷卻至室溫時包殼直徑軸向分布情況。由于BISON程序采用MATPRO塑性模型[24]，而BEEs和FRAPTRAN程序均采用保守的PNNL塑性模型，因此兩者計算的包殼變形高于BISON程序。此外，由于BEEs程序考慮了芯塊的彈塑性變形和蠕變行為而FRAPTRAN程序模擬時將芯塊作為剛體，導(dǎo)致BEEs程序計算值略低于FRAPTRAN程序。

表8列出Na-2實驗各程序計算值與實驗值對比。由表8和圖11可知，各程序計算值均總體低于實驗測量值，這主要是由于缺乏冷卻劑及包殼約束的相關(guān)實驗數(shù)據(jù)，各程序模擬時邊界條件設(shè)置有偏差。此外，各程序計算PCMI時采用的無摩擦接觸模型忽略了應(yīng)力的軸向分量，使包殼有效應(yīng)力計算偏小，進一步導(dǎo)致包殼變形預(yù)測值偏低。芯塊兩端的碟形倒角使包殼在PCMI中產(chǎn)生較大的局部變形，導(dǎo)致實驗測量值在軸向會呈現(xiàn)局部波動，而考慮到缺乏實驗芯塊倒角的詳細描述和PCMI計算的復(fù)雜性，當(dāng)前各程序模擬過程中并未對碟形倒角進行建模，后續(xù)可通過更精確的二維、三維建模對該現(xiàn)象進行分析。

表8 Na-2實驗各程序計算值與實驗值對比Table 8 Comparison among code calculation value and test data for Na-2

圖11 包殼直徑對比Fig.11 Comparison of cladding diameter

3.2 LOCA模擬驗證

NRU-LOCA MT4實驗[25]是PNNL在加拿大Chalk River國家實驗室的NRU反應(yīng)堆中開展的LOCA瞬態(tài)實驗。實驗中將零功率燃料棒進行輻照、升溫、再淹沒，測量了包殼內(nèi)表面溫度，實驗后拆除實驗段，測量了燃料棒變形。實驗中燃料棒設(shè)計參數(shù)列于表9，燃料類型為UO2，包殼類型為Zr4。本文將BEEs程序與FRAPTRAN程序[26]和BISON程序[17]計算結(jié)果進行對比。模擬時網(wǎng)格劃分如圖12所示，采用網(wǎng)格類型為QUAD4，在燃料徑向劃分11個節(jié)點，軸向劃分1 680個節(jié)點，包殼徑向劃分4個節(jié)點，軸向劃分2 500個節(jié)點。

圖12 MT4實驗算例網(wǎng)格劃分Fig.12 Mesh partition of MT4 test case

表9 MT4實驗設(shè)計參數(shù)Table 9 Design parameter of MT4 test

圖13示出中心高度處(183 cm)包殼內(nèi)表面溫度的對比。前10 s實驗維持穩(wěn)態(tài)工況，之后發(fā)生LOCA使溫度上升。值得注意的是，由于FRAPTRAN程序軸向節(jié)點劃分較少，其計算值取自離軸向高度183 cm最近的節(jié)點。由圖13可知，BEEs程序計算結(jié)果與實驗值符合良好，與FRAPTRAN和BISON程序的結(jié)果基本一致。

圖13 包殼內(nèi)表面中心高度處溫度的變化Fig.13 Change of temperature of cladding inner surface at axial mid-plane

燃料棒內(nèi)壓變化對比如圖14所示。由于BEEs程序采用與FRAPTRAN程序一致的溫度邊界條件，使燃料棒溫度計算偏高，進而導(dǎo)致內(nèi)壓計算值高于實驗值。針對LOCA下包殼高溫變形的模擬，BEEs和BISON程序均采用連續(xù)蠕變行為模型，因此兩者計算的內(nèi)壓變化趨勢總體一致。而FRAPTRAN程序需在基于小變形的彈塑性求解及單獨的鼓脹變形模塊BALON2間進行模型轉(zhuǎn)換[27]，導(dǎo)致其內(nèi)壓計算曲線中出現(xiàn)明顯轉(zhuǎn)折點和兩個峰值。

圖14 燃料棒內(nèi)壓的變化Fig.14 Change of rod internal pressure

圖15為失效時燃料溫度分布云圖。由于LOCA條件下的高溫蠕變效應(yīng)，在失效時刻包殼中間高度處會發(fā)生顯著的鼓脹現(xiàn)象。

包殼破裂失效時的參數(shù)對比列于表10。對于包殼失效時間、失效內(nèi)壓和失效溫度，BEEs與BISON程序計算值及實驗值符合良好。BEEs與FRAPTRAN程序的誤差主要由高溫下包殼變形的模擬方法不同所導(dǎo)致，F(xiàn)RAPTRAN程序預(yù)測的失效參數(shù)過于保守。由表10可知，各程序計算的破裂時刻包殼環(huán)向應(yīng)變相當(dāng)，但由于上述程序目前均無法模擬破裂后包殼破口的真實形態(tài)，而實驗數(shù)據(jù)是對包殼破口進行實際測量所得，因此模擬結(jié)果與實驗數(shù)據(jù)有較大偏差。

圖中將軸向高度縮短為原長0.2倍，徑向長度增大為原長15倍圖15 包殼失效時燃料溫度分布Fig.15 Distribution of fuel temperature when cladding burst

表10 包殼失效參數(shù)對比Table 10 Comparison of parameter when cladding burst

4 總結(jié)

本文對基于MOOSE平臺的棒狀燃料元件性能分析程序BEEs進行瞬態(tài)功能開發(fā)，并對開發(fā)的RIA和LOCA模塊開展模塊驗證及整體驗證。結(jié)果表明，BEEs程序能實現(xiàn)瞬態(tài)事故條件下燃料溫度、力學(xué)變形、內(nèi)壓演變的模擬及包殼失效的預(yù)測，計算結(jié)果合理，并得到以下結(jié)論。

1) 針對瞬態(tài)功能擴展模型的驗證表明：BEEs程序中添加的塑性求解模塊在均勻塑性變形和強化階段相比于實驗值的最大相對誤差分別為4.7%和11.3%；高溫蠕變模型適用于環(huán)向應(yīng)變大于1%的情況，最大相對誤差為8.6%；在0.1～100 K·s-1升溫速率范圍內(nèi)，包殼相變模型計算值的相對誤差均在1.85%以內(nèi)；高溫氧化模型適用于溫度較高(大于973 K)、氧化層較厚(大于10 μm)的情況，最大相對誤差為11.4%。

2) 針對CABRI-RIA Na-2實驗，BEEs程序計算值與其他程序符合良好，與FRAPTRAN程序相比，芯塊和包殼溫度計算最大相對誤差分別為0.4%和2.0%，包殼環(huán)向應(yīng)變和殘余變形相對誤差分別為1.2%和1.8%。由于缺乏冷卻劑及包殼約束的相關(guān)實驗數(shù)據(jù)和無摩擦接觸模型的局限性，導(dǎo)致模擬PCMI時程序計算的殘余變形與實驗值有較大偏差。

3) 針對NRU-LOCA MT4實驗，BEEs程序計算值模擬結(jié)果具有合理性，與BISON和FRAPTRAN程序?qū)Ρ惹闆r良好，BEEs程序預(yù)測的包殼失效時間和內(nèi)壓均在實驗值范圍內(nèi)，失效溫度計算值相對誤差在0.4%～3.7%之間。由于程序目前無法模擬到達破裂限值時包殼實際的破口形態(tài)，因此針對破口處環(huán)向應(yīng)變的計算偏差較大。