邵世豪,劉宙宇,許曉北,宗育凡,曹良志,吳宏春

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院,陜西 西安 710049)

較高的鈾資源利用率和較硬的中子能譜等優(yōu)點(diǎn)使快堆技術(shù)成為我國核電“三步走”戰(zhàn)略中的重要一步。液態(tài)金屬冷卻快堆作為第四代核能系統(tǒng)極具潛力的堆型之一,對(duì)于我國核能創(chuàng)新發(fā)展具有非常重要的意義[1]。

燃料元件是核反應(yīng)堆最基本的部件之一,對(duì)其進(jìn)行的燃料元件性能分析是反應(yīng)堆研發(fā)中最具挑戰(zhàn)性的工作之一。燃料元件性能分析用于預(yù)測燃料行為與關(guān)鍵參數(shù)的演化,是指導(dǎo)燃料元件設(shè)計(jì)、預(yù)測燃料元件服役壽命、確保反應(yīng)堆安全運(yùn)行的重要課題之一。國內(nèi)針對(duì)壓水堆燃料性能分析程序已經(jīng)開展了大量研究,然而,相較于壓水堆,快堆的運(yùn)行環(huán)境更為嚴(yán)苛,其通量水平更高,運(yùn)行溫度也更高,這導(dǎo)致堆芯材料結(jié)構(gòu)在快堆中受到的輻照影響比壓水堆更為顯著[2]。因此,開發(fā)一套適用于快堆環(huán)境的燃料性能分析程序,對(duì)于優(yōu)化快堆設(shè)計(jì)、提升其安全性以及確保反應(yīng)堆的長期穩(wěn)定運(yùn)行具有至關(guān)重要的意義。

從20世紀(jì)開始,許多國家相繼開發(fā)了一系列針對(duì)液態(tài)金屬冷卻快堆燃料元件性能分析程序,如美國麻省理工學(xué)院(MIT)開發(fā)的FEAST-OXIDE[3],適用于氧化物燃料和金屬燃料,被認(rèn)為是目前考慮物理現(xiàn)象最為全面的快堆燃料性能分析程序;法國原子能與可替代能源委員會(huì)(CEA)開發(fā)的GERMINAL[4]以及日本原子能委員會(huì)開發(fā)的FEMAXI-FBR[5]等。國內(nèi)近些年也開發(fā)了一些快堆燃料元件性能分析軟件,如中國原子能科學(xué)研究院開發(fā)的LIFEANLS[6],是國內(nèi)第1個(gè)快堆燃料性能分析程序,以及后續(xù)開發(fā)的FIBER程序[7-8];中國科學(xué)技術(shù)大學(xué)開發(fā)的KMC-FUEL[2]等。

為實(shí)現(xiàn)快堆燃料元件性能分析程序自主化,本文基于多物理耦合平臺(tái)MOOSE,開發(fā)燃料元件性能分析程序LoongCALF,該程序旨在模擬燃料元件在穩(wěn)態(tài)工況下的長壽期服役性能演化,為LMFR燃料元件設(shè)計(jì)、服役壽命預(yù)測及安全評(píng)價(jià)提供科學(xué)依據(jù)。本文首先對(duì)LoongCALF程序進(jìn)行介紹,隨后,通過構(gòu)建并計(jì)算兩個(gè)設(shè)計(jì)算例,將LoongCALF程序與中國原子能科學(xué)研究院開發(fā)的Fiber-Oxide程序進(jìn)行深入對(duì)比分析,確保LoongCALF程序的可靠性和有效性,為其在實(shí)際應(yīng)用中的推廣使用提供堅(jiān)實(shí)的基礎(chǔ)。

1 方法

1.1 程序總體介紹

LoongCALF程序是在之前的輕水堆燃料性能分析程序NECP-CALF[9]的基礎(chǔ)上,擴(kuò)展而來的針對(duì)液態(tài)金屬冷卻快堆燃料元件的燃料性能分析程序。程序的體系結(jié)構(gòu)分為4個(gè)層次:控制層、求解層、物理層與材料層,如圖1所示。其中:控制層負(fù)責(zé)軟件的啟動(dòng)、執(zhí)行與調(diào)度;求解層負(fù)責(zé)熱-力耦合方程的求解;物理層負(fù)責(zé)燃料中物理現(xiàn)象及過程的模擬;材料層負(fù)責(zé)計(jì)算材料行為參數(shù)。

圖1 程序體系結(jié)構(gòu)示意圖Fig.1 Schematic diagram of program architecture

程序的主要處理流程如圖2所示。

圖2 程序主要處理流程圖Fig.2 Main processing flowchart of program

在燃料性能分析中,熱-力耦合問題涉及溫度場和應(yīng)力場之間的動(dòng)態(tài)交互,這些交互通過材料的熱膨脹系數(shù)、熱傳導(dǎo)率等物理參數(shù)得以體現(xiàn)。針對(duì)這一問題,程序采用Picard迭代方法(又稱固定點(diǎn)迭代方法)進(jìn)行求解。該方法將多物理場問題分解為多個(gè)獨(dú)立的子問題,這些子問題(如溫度場、應(yīng)力場等)通過MOOSE框架中的MultiApp模塊進(jìn)行模擬求解。在求解某個(gè)物理場過程中,其他物理場采用初始值或者上一迭代步的計(jì)算結(jié)果,通過迭代更新耦合物理變量,直至滿足預(yù)設(shè)的收斂條件。Picard迭代方法的應(yīng)用顯著提升了燃料在復(fù)雜熱力環(huán)境下的性能預(yù)測的準(zhǔn)確性和效率。

對(duì)于單獨(dú)的熱學(xué)或力學(xué)問題求解,JFNK方法(Jacobian-free Newton-Krylov method)已成為一種高效且實(shí)用的數(shù)值求解策略。該方法巧妙地結(jié)合了Newton法的收斂性與Krylov子空間方法的計(jì)算效率,同時(shí)避免了顯式計(jì)算或存儲(chǔ)完整雅可比矩陣的需求,從而顯著降低了計(jì)算資源和內(nèi)存的使用量。然而,為進(jìn)一步提升JFNK方法的性能和穩(wěn)定性,通常采用預(yù)處理計(jì)算方法來優(yōu)化其求解過程。預(yù)處理矩陣的引入旨在改善線性系統(tǒng)的條件數(shù),進(jìn)而加快迭代收斂速度。在這一過程中,基于MOOSE框架本身的預(yù)處理功能,而MOOSE則采用解析表達(dá)式來定義預(yù)處理矩陣。這種方法不僅提高了求解的精度,還增強(qiáng)了算法的魯棒性,使得能夠更加準(zhǔn)確地模擬和預(yù)測物理現(xiàn)象。

1.2 基于MOOSE平臺(tái)程序開發(fā)

MOOSE(multiphysics object oriented simulation environment)平臺(tái)是一個(gè)開源的,能夠同時(shí)模擬多個(gè)物理場,用于構(gòu)建復(fù)雜系統(tǒng)的仿真模型。它支持多維度建模計(jì)算,可同時(shí)模擬多個(gè)物理場,如熱傳導(dǎo)與力學(xué)物理場耦合計(jì)算。在MOOSE中,可使用模塊化的方式來創(chuàng)建復(fù)雜的仿真模型,每個(gè)模塊可代表1個(gè)物理場或1個(gè)子系統(tǒng),通過組合這些模塊,可構(gòu)建出多維度、多物理場的仿真模型。此外,MOOSE還提供了豐富的物理模型庫,包括常見的材料模型、邊界條件、初始條件等,方便用戶快速構(gòu)建仿真模型。

表1 燃料元件結(jié)構(gòu)參數(shù)Table 1 Structural parameter of fuel element

程序基于MOOSE平臺(tái)進(jìn)行開發(fā),充分利用了MOOSE平臺(tái)的優(yōu)勢。程序采用模塊化設(shè)計(jì),有單獨(dú)的材料模塊,目前燃料類型支持UO2及MOX,包殼材料支持HT9及1515Ti,冷卻劑材料支持鈉、鉛及鉛鉍。此外程序還具備向316(Ti)SS、T92等先進(jìn)包殼材料及U-Zr、U-Pu-Zr金屬燃料擴(kuò)展的架構(gòu)條件。

為滿足不同復(fù)雜度和需求的模型計(jì)算,程序?yàn)橛脩籼峁┝硕喾N建模方式和多維度網(wǎng)格求解功能。建模方式包括可視化建模和基于腳本的數(shù)字化建模兩種方式?？梢暬７绞街庇^易用,用戶可通過簡單的拖拽和配置即可完成模型的構(gòu)建;而基于腳本的數(shù)字化建模方式則更加靈活,用戶可通過編寫腳本來自定義模型和計(jì)算流程。此外,程序支持一維、二維和三維網(wǎng)格建模計(jì)算,以滿足不同尺寸和形狀的燃料元件模型構(gòu)建需求。目前,對(duì)于燃料元件建模,程序通常采用1.5維建模方法。這種方式不僅大幅減少了計(jì)算網(wǎng)格的數(shù)量和復(fù)雜度,簡化了模型構(gòu)建和維護(hù)過程,還降低了對(duì)計(jì)算資源和時(shí)間的需求。同時(shí),它依然能夠捕捉燃料行為的關(guān)鍵特征,為分析提供有效的支持。

1.3 數(shù)學(xué)模型

1) 控制方程

燃料元件的溫度通過熱傳導(dǎo)方程求解:

(1)

式中:ρ為密度,kg/m3;cp為比熱,J/(kg·℃);T為溫度,K;k為熱導(dǎo)率;S為熱源。

燃料元件的應(yīng)力-應(yīng)變狀態(tài)通過固體力學(xué)平衡方程求解:

(2)

式中:σ為應(yīng)力;f為單位質(zhì)量的體積力。

幾何方程為:

(3)

式中:ε為應(yīng)變;u為位移,m。

本構(gòu)方程為:

σ=E:(ε-εplastic-εcreep-εeigen)

(4)

式中:E為彈性張量;εplastic為塑性應(yīng)變;εcreep為蠕變應(yīng)變;εeigen為本征應(yīng)變。

2) 流體換熱模型

程序采用單通道假設(shè)[2],燃料元件與環(huán)繞它的流體通道換熱,通道與通道之間沒有物質(zhì)和能量的交流。冷卻劑溫度控制方程為:

(5)

式中:c為冷卻劑比熱容,J/(kg·K);T為冷卻劑溫度,K;w為冷卻劑軸向流速,m/s;qv,ch為冷卻劑流道中的體積熱源,W/m3,包括γ輻射帶出的少量沉積在冷卻劑區(qū)域的裂變能及包殼向冷卻劑傳遞的熱量(以體積熱源的形式引入)。

3) 裂變氣體行為模型

模擬燃料中裂變氣體的釋放常用的方法為兩步法:第1步,裂變氣體在晶粒內(nèi)部產(chǎn)生并擴(kuò)散到晶界,在晶界處逐漸形成氣泡;第2步,裂變氣體的氣泡在晶界處不斷長大、合并,最終與外界氣體環(huán)境連通,使得儲(chǔ)存在晶界中的裂變氣體被釋放到外界氣體環(huán)境中[10]。

第1階段的氣體擴(kuò)散過程由下式描述[11]:

(6)

式中:Cg為晶內(nèi)裂變氣體濃度;t為時(shí)間;Deff為有效氣體擴(kuò)散系數(shù);r為一維球坐標(biāo)系的坐標(biāo);β為氣體產(chǎn)生速率。

第2階段,當(dāng)晶界氣體與外界連通(即晶界氣體達(dá)到飽和)時(shí),其在晶界上的面濃度由下式計(jì)算:

(7)

f(θdihe)=1-1.5cosθdihe+0.5cos3θdihe

(8)

式中:Ns為裂變氣體在晶界上的飽和面濃度,m-2;rb為晶界處裂變氣體氣泡半徑,m;θdihe為晶界氣泡的二面角;Vc為飽和時(shí)晶界被氣泡覆蓋的比例;kB為玻耳茲曼常數(shù),kB=1.380 6×10-23J/K;γ為氣泡表面張力,J/m2;pext為外部靜水壓力,Pa。

晶界處的裂變氣體積攢超過限值時(shí),多余的氣體將被釋放到燃料元件氣腔中。

4) 燃料化學(xué)與結(jié)構(gòu)調(diào)整

(1) 氧金屬比-燃耗依賴關(guān)系

燃料中的氧金屬比會(huì)影響燃耗的速度和程度,通過研究氧金屬比與燃耗的相關(guān)性,可更準(zhǔn)確地預(yù)測燃料的燃燒壽命和性能。

氧金屬比的計(jì)算參考文獻(xiàn)[12],即:

(9)

對(duì)于VU=4,且VPu<4,有:

(2YBa-Sr+4YZr-Nb+3YY-Re+4fMoYMo)β

(10)

對(duì)于VU>4,且VPu=4,有:

(2YBa-Sr+4YZr-Nb+3YY-Re+4fMoYMo)β

(11)

控制燃料的初始氧金屬比,使鈾的價(jià)態(tài)在整個(gè)輻照循環(huán)中不會(huì)顯著超過4,這將限制燃料包殼化學(xué)相互作用(FCCI)。

(2) 孔隙遷移

通常,氧化物燃料中包含了制造過程中形成的孔隙,這些孔隙被氦氣所填充。由于這些孔隙體積大且保持低壓,因此其形狀呈透鏡狀。當(dāng)燃料處于運(yùn)行功率下,這些孔隙會(huì)沿著溫度梯度向燃料的內(nèi)部區(qū)域遷移。這一遷移機(jī)制主要基于蒸汽傳輸過程:孔隙在熱表面發(fā)生蒸發(fā),并在冷表面發(fā)生冷凝。因此,孔隙會(huì)逐漸遷移到燃料棒的中心區(qū)域,最終形成一個(gè)中心空洞[3]。

孔隙遷移速度按照以下關(guān)系式計(jì)算:

Vp=0.337 6exp(-66 490/T)

(12)

式中,Vp為孔隙速度,m/s。

孔隙率分布函數(shù)Np(r,t)是用于描述孔隙率隨時(shí)間以及燃料棒中徑向位置而變化的數(shù)學(xué)表達(dá)式,具體可由以下方程求解:

(13)

Np(r,t=0)=Np0

(14)

Np(r=R,t)=Np0

(15)

式中:Np0為制造的燃料中的孔隙濃度;Np(r,t=0)描述方程求解的初值;Np(r=R,t)描述方程求解的邊界條件。

(3) 晶粒生長模型

在反應(yīng)堆中,燃料的晶粒會(huì)逐漸生長,晶粒的尺寸是裂變氣體行為模型中的一個(gè)重要參數(shù),晶粒生長模型采用Ainscough等[13]提出,Khoruzhii等[14]修正的晶粒長大動(dòng)力學(xué)[2],即:

(16)

式中:a為晶粒直徑,μm;K為晶粒生長率,μm2/h;amax為可穩(wěn)定存在的晶粒上限尺寸,μm;airr為考慮輻照效應(yīng)的修正,μm。

5) 材料行為與物性模型

(1) 熱導(dǎo)率

MOX燃料熱導(dǎo)率由下式計(jì)算[3]:

k=F1F2F3F4ko

(17)

(18)

(19)

(20)

(21)

F4=1-αP

(22)

式中:k為MOX燃料熱導(dǎo)率;β為燃耗,at%;P為燃料孔隙率;α為孔隙系數(shù)(1.5或2.5),本文取α=2.5。

(2) 熱膨脹

熱膨脹模型對(duì)燃料芯塊的尺寸變化進(jìn)行建模。燃料的尺寸變化會(huì)影響芯塊到包殼的間隙尺寸,這是決定間隙傳熱的一個(gè)主要因素,從而決定儲(chǔ)存能量,這是安全分析的一個(gè)重要量。MOX燃料的線熱膨脹應(yīng)變模型[3]為:

(23)

式中:ED為缺陷生成能,J;k=1.38×10-23J/K;K1,2,3為模型系數(shù)。

(3) 重定位

當(dāng)核反應(yīng)堆運(yùn)行時(shí),芯塊會(huì)受到高溫和高壓的作用,導(dǎo)致燃料碎片的熱膨脹和震動(dòng)移位。MOX燃料直徑遷移由下式計(jì)算:

DR=0.111GDci-45

(24)

式中:DR為直徑遷移,μm;G為直徑間隙,μm;Dci為包殼內(nèi)徑,mm。

為保證芯塊外徑與重定位模型計(jì)算值一致,令芯塊中央空洞也產(chǎn)生重定位應(yīng)變。重定位應(yīng)變只產(chǎn)生在R-θ平面上,其線應(yīng)變?yōu)?

εl=DR/2Rpout

(25)

(4) 腫脹

MOX燃料裂變反應(yīng)產(chǎn)生的裂變產(chǎn)物大部分為固態(tài),小部分為氣態(tài),這些產(chǎn)物滯留在燃料中,造成燃料的腫脹。

固態(tài)裂變產(chǎn)物導(dǎo)致的腫脹被稱為固態(tài)腫脹,可由下式[15]計(jì)算:

Δεsw-s=5.577×10-5ρΔBu

(26)

式中:Δεsw-s為固態(tài)腫脹體應(yīng)變?cè)隽?ρ為密度,kg/m3;ΔBu為燃耗增量。

氣態(tài)裂變產(chǎn)物導(dǎo)致的腫脹被稱為氣態(tài)腫脹,可由下式計(jì)算:

Δεsw-g=1.96×10-31ρΔBu(2 800-T)11.73×

exp(-0.016 2(2 800-T)-0.017 8ρBu)

(27)

式中,Δεsw-g為氣態(tài)腫脹體應(yīng)變?cè)隽俊?/p>

(5) 密實(shí)化

燃料芯塊在服役過程中,隨著裂變產(chǎn)物的累積和溫度的提升,燃料中的孔隙不斷收縮,導(dǎo)致芯塊整體的收縮,這樣的現(xiàn)象被稱為密實(shí)化。

MOX燃料密度分?jǐn)?shù)[3]由下式計(jì)算:

ρ=ρo+(0.96-ρo)(1-exp(-Bu/0.6))

(28)

式中:ρ為燃料密度分?jǐn)?shù);ρo為初始燃料密度分?jǐn)?shù)。

密實(shí)化導(dǎo)致的應(yīng)變:

εV=ρo/ρ-1

(29)

(6) 蠕變

燃料芯塊及包殼在高溫情況下,即使其應(yīng)力小于材料的屈服極限,也會(huì)緩慢地發(fā)生不可逆的形變,這種現(xiàn)象叫做蠕變。

MOX燃料熱蠕變由下式計(jì)算[3]:

(30)

輻照蠕變:

(31)

2 程序驗(yàn)證結(jié)果

為驗(yàn)證本文開發(fā)的液態(tài)金屬冷卻快堆燃料元件性能分析程序LoongCALF的計(jì)算準(zhǔn)確性及適用性,設(shè)計(jì)了兩個(gè)基準(zhǔn)算例進(jìn)行數(shù)值模擬與計(jì)算。算例1是以UO2為芯塊、1515Ti為包殼的燃料棒;算例2是以MOX為芯塊、HT9為包殼的燃料棒。Fiber-Oxide程序是由中國原子能科學(xué)研究院開發(fā)的鈉冷快堆燃料元件性能分析程序[7],被選為本文的參考程序。兩個(gè)算例的燃料元件結(jié)構(gòu)參數(shù)列于表1,沿軸向功率分布如圖3所示。

圖3 軸向功率分布Fig.3 Axial power distribution

使用LoongCALF程序?qū)蓚€(gè)基準(zhǔn)算例進(jìn)行數(shù)值模擬與計(jì)算,將LoongCALF程序計(jì)算得到的芯塊最高溫度、間隙寬度、應(yīng)力、氣腔壓強(qiáng)及裂變氣體釋放相關(guān)參數(shù)與Fiber-Oxide程序計(jì)算結(jié)果進(jìn)行比較,本文僅展示最大功率的軸向段5的各參數(shù)隨時(shí)間點(diǎn)變化,其他功率軸向段結(jié)論相同。

2.1 芯塊最高溫度對(duì)比分析

兩個(gè)算例軸向段5芯塊最高溫度比較結(jié)果如圖4所示。

圖4 燃料最高溫度比較Fig.4 Comparison of the highest fuel temperature

燃料性能溫度分析模塊中,芯塊包殼間隙寬度對(duì)溫度產(chǎn)生直接且重大的影響。這種影響主要體現(xiàn)在間隙寬度存在差異的情況下,芯塊內(nèi)表面溫度會(huì)呈現(xiàn)出顯著的差異。此外,芯塊熱導(dǎo)率模型的差異也是導(dǎo)致溫度差異的重要因素。參考文獻(xiàn)[16]中各國際知名燃料性能分析程序計(jì)算得到的溫度差異最大在300 K左右,對(duì)于兩個(gè)算例,兩個(gè)程序的最大絕對(duì)偏差為150 K左右,最大相對(duì)偏差為10%以內(nèi),在可接受范圍之內(nèi)。

2.2 間隙寬度對(duì)比分析

兩個(gè)算例軸向段5間隙寬度比較結(jié)果如圖5所示。

圖5 間隙寬度比較Fig.5 Comparison of gap width

間隙寬度的初始差異源于兩個(gè)程序在處理重定位應(yīng)變時(shí)的理念不同。兩個(gè)程序?qū)τ谛緣K包殼接觸時(shí)間預(yù)測大致相同,但由于影響芯塊包殼間隙的物理量諸多,如熱膨脹、腫脹、密實(shí)化、蠕變等都會(huì)對(duì)芯塊包殼作用產(chǎn)生應(yīng)變,不同程序間很難做到完全一致,總體符合較好。

2.3 應(yīng)力對(duì)比分析

兩個(gè)算例軸向段5包殼內(nèi)表面米塞斯應(yīng)力比較結(jié)果如圖6所示。

圖6 包殼內(nèi)表面米塞斯應(yīng)力比較Fig.6 Comparison of cladding inner surface Mises stress

包殼應(yīng)力的不同主要源于兩種方法對(duì)包殼處理方式的差異。Fiber-Oxide采用了廣義平面應(yīng)變假設(shè),而LoongCALF則采用二維計(jì)算方式。但可看到兩個(gè)程序?qū)?yīng)力趨勢的預(yù)測大致相同,符合較好。

2.4 壓強(qiáng)及裂變氣體對(duì)比分析

兩個(gè)算例氣腔壓強(qiáng)及裂變氣體釋放份額比較結(jié)果如圖7、8所示。

圖7 氣腔壓強(qiáng)比較Fig.7 Comparison of plenum pressure

圖8 裂變氣體釋放份額比較Fig.8 Comparison of fission gas release fraction

氣腔壓強(qiáng)受裂變氣體影響,由于兩個(gè)程序裂變氣體開始釋放的時(shí)間不同而略有差距。Fiber-Oxide程序計(jì)算得到的裂變氣體開始釋放的燃耗明顯早于LoongCALF程序,但考慮到裂變氣體釋放計(jì)算具有較大不確定度以及不同模型之間的差異,這樣的結(jié)果是可接受的,同時(shí)觀察到,兩個(gè)程序?qū)\(yùn)行后期的裂變氣體釋放份額模擬基本相同,符合較好。

3 結(jié)論

本文開發(fā)的快堆燃料元件性能分析程序LoongCALF具備對(duì)快堆UO2及MOX燃料的一系列分析功能,包括導(dǎo)熱、間隙傳熱、重結(jié)構(gòu)、裂變產(chǎn)物遷移、裂變氣體釋放、輻照腫脹、輻照蠕變等計(jì)算功能。為驗(yàn)證程序的準(zhǔn)確性,本文設(shè)計(jì)了兩個(gè)完整燃料棒算例進(jìn)行程序?qū)Ρ葴y試。結(jié)果表明,LoongCALF程序能夠準(zhǔn)確模擬液態(tài)金屬冷卻快堆穩(wěn)態(tài)工況條件下燃料元件內(nèi)部的燃料行為與關(guān)鍵參數(shù)演化。

鑒于燃料行為涉及多種復(fù)雜現(xiàn)象的耦合,且不同模型間的兼容性對(duì)分析結(jié)果至關(guān)重要,未來將依托更多實(shí)際測量案例,對(duì)LoongCALF程序進(jìn)行更深入的驗(yàn)證,以確保其準(zhǔn)確性和合理性。同時(shí),程序可基于2.5維和三維建模開展更多的復(fù)雜幾何工況分析。

本研究得到了中國原子能科學(xué)研究院Fiber-Oxide程序結(jié)果的幫助,為本研究提供了重要的比較數(shù)據(jù)。在此,對(duì)中國原子能科學(xué)研究院的貢獻(xiàn)表示衷心的感謝,同時(shí)感謝陳啟董老師的經(jīng)驗(yàn)分享與技術(shù)指導(dǎo)。