基于MOOSE平臺的棒狀燃料元件性能瞬態分析程序開發與驗證

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

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

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

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

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

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

1 數學物理模型

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

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

1.1 彈塑性模型[13]

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

σ=ε·E

(1)

其中：σ為應力，MPa；E為彈性模量，MPa；ε為應變。

當包殼的應力達到屈服極限后，用冪次定律描述包殼的塑性變形：

(2)

1.2 高溫蠕變模型

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

(3)

表1 Zr4包殼蠕變常數Table 1 Creep constant of Zr4 cladding

1.3 高溫相變模型

隨著包殼溫度上升，鋯合金晶體從六角α相結構轉變為立方形β相結構。鋯包殼的相變速率可用下式[15]表示：

(4)

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

(5)

(6)

其中：w為局部氫濃度，ppm；k0為動力常數，s-1；E為等效激活能，J；kb為玻爾茲曼常數，J·K-1；km為常數，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。

常數km取值為：

(7)

相變轉換起始溫度為：

(8)

(9)

1.4 高溫氧化模型

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

(10)

(11)

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

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

表2 氧化關系式參數Table 2 Parameter of oxidation relationship

1.5 包殼失效準則[14]

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

1) 包殼局部應力限制準則

(12)

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

表3 局部應力限制準則參數Table 3 Parameter of local stress limiting criterion

2) 塑性“失穩”準則

(13)

2 瞬態模型驗證

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

2.1 包殼彈塑性求解驗證

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

表4 應力應變關系式驗證算例參數Table 4 Parameter of validation case of stress-strain correlation

2.2 高溫蠕變模型驗證

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

圖3 內壓1.4 MPa時包殼環向應變對比Fig.3 Comparison of cladding strain with internal pressure at 1.4 MPa

表5 Hardy tube實驗設計參數Table 5 Design parameter of Hardy tube test

2.3 包殼相變模型驗證

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

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

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

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

3 整體驗證分析

3.1 RIA模擬驗證

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

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

圖8 Na-2實驗算例網格劃分Fig.8 Mesh partition of case Na-2

表7 Na-2實驗設計參數Table 7 Design parameter of Na-2 test

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

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

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

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

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

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

表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反應堆中開展的LOCA瞬態實驗。實驗中將零功率燃料棒進行輻照、升溫、再淹沒，測量了包殼內表面溫度，實驗后拆除實驗段，測量了燃料棒變形。實驗中燃料棒設計參數列于表9，燃料類型為UO2，包殼類型為Zr4。本文將BEEs程序與FRAPTRAN程序[26]和BISON程序[17]計算結果進行對比。模擬時網格劃分如圖12所示，采用網格類型為QUAD4，在燃料徑向劃分11個節點，軸向劃分1 680個節點，包殼徑向劃分4個節點，軸向劃分2 500個節點。

圖12 MT4實驗算例網格劃分Fig.12 Mesh partition of MT4 test case

表9 MT4實驗設計參數Table 9 Design parameter of MT4 test

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

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

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

圖14 燃料棒內壓的變化Fig.14 Change of rod internal pressure

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

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

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

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

4 總結

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

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

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

3) 針對NRU-LOCA MT4實驗，BEEs程序計算值模擬結果具有合理性，與BISON和FRAPTRAN程序對比情況良好，BEEs程序預測的包殼失效時間和內壓均在實驗值范圍內，失效溫度計算值相對誤差在0.4%～3.7%之間。由于程序目前無法模擬到達破裂限值時包殼實際的破口形態，因此針對破口處環向應變的計算偏差較大。