壓水堆燃料輻照性能多物理耦合并行分析程序開發(fā)

韓智杰，何曉軍,*，明 春，楊衍康，任 帥，胡長軍，楊 文

(1.中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究所，北京 102413；2.北京科技大學(xué)，北京 100083)

核燃料元件是反應(yīng)堆系統(tǒng)中的核心部件，是反應(yīng)堆運(yùn)行過程中熱量及裂變反應(yīng)的主要來源，故燃料元件的堆內(nèi)性能對反應(yīng)堆經(jīng)濟(jì)性和安全性有重要影響。正常運(yùn)行條件下，燃料元件設(shè)計(jì)必須保證其在高溫、高壓、強(qiáng)輻照環(huán)境下的完整性。利用計(jì)算模擬的方法可以用較低的成本預(yù)測燃料元件堆內(nèi)行為，達(dá)到評(píng)價(jià)元件堆內(nèi)輻照服役性能、指導(dǎo)燃料元件設(shè)計(jì)的目的[1]。因此，燃料元件性能分析程序?qū)θ剂涎邪l(fā)[2]、安全評(píng)價(jià)具有重要意義[3]。

燃料元件的堆內(nèi)行為涉及中子物理、熱工水力、固體力學(xué)、化學(xué)腐蝕等多種物理現(xiàn)象[4]，典型商用壓水堆中裝載上萬根燃料元件，傳統(tǒng)分析評(píng)價(jià)方法通常選取典型或極限工況燃料元件進(jìn)行分析，存在一定的局限性，故有必要建立先進(jìn)多物理耦合并行化燃料性能分析軟件進(jìn)行多棒并行分析，高效率、高精度評(píng)價(jià)反應(yīng)堆燃料性能，為提高反應(yīng)堆安全性及經(jīng)濟(jì)性提供有效依據(jù)[5]。

本工作建立燃料元件溫度、變形、裂變氣體釋放及內(nèi)壓等主要分析模型，通過熱工-力學(xué)-內(nèi)壓多物理耦合，基于先進(jìn)并行計(jì)算方法，開發(fā)先進(jìn)多物理燃料性能分析程序Athena。

1 數(shù)理模型

1.1 程序介紹

Athena是中國數(shù)值反應(yīng)堆原型系統(tǒng)(CVR1.0)的重要組成程序，CVR1.0目標(biāo)是開發(fā)面向E級(jí)超算的反應(yīng)堆模擬軟件，包括中子物理、熱工水力、結(jié)構(gòu)力學(xué)、燃料性能和材料性能等核心分析軟件及多物理耦合環(huán)境[6]。為充分發(fā)揮超算優(yōu)勢，適應(yīng)全堆芯pin-by-pin[7]物理-熱工-燃料精細(xì)化多物理耦合應(yīng)用需求[8]，Athena利用先進(jìn)并行算法，建立了燃料并行分析能力，具有高效、精細(xì)化全堆芯燃料性能分析優(yōu)勢。

1.2 溫度分析模型

燃料元件溫度計(jì)算在性能分析中占有重要地位，它是所有分析的基礎(chǔ)。溫度分析模型的主要功能是獲得燃料元件的溫度分布，即計(jì)算熱量在冷卻劑、包殼、芯塊-包殼間隙、芯塊內(nèi)部的傳遞。

將燃料元件及冷卻劑簡化為單通道單相傳熱條件，冷卻劑溫度通過求解質(zhì)量、動(dòng)量和能量守恒方程得到。冷卻劑到包殼傳熱方式為對流換熱，包殼表面溫度可通過牛頓冷卻定律得到[9]。

包殼及燃料芯塊可簡化為不同熱導(dǎo)率材料組成的軸對稱圓柱體，故其內(nèi)部熱量傳遞可用穩(wěn)態(tài)傅里葉導(dǎo)熱方程(式(1))描述，差分求解節(jié)點(diǎn)劃分示意圖如圖1所示。圖1中：δ為徑向節(jié)點(diǎn)間距，m；N為徑向節(jié)點(diǎn)數(shù)量。

圖1 燃料棒溫度差分求解節(jié)點(diǎn)示意圖Fig.1 Schematic of fuel rod node for temperature model

(1)

式中：k為熱導(dǎo)率，W/(m·K)；S為體積元表面積，m2；n為表面標(biāo)準(zhǔn)單位向量；q(r)為熱源，W/m3；T為溫度，K；V為體積元體積，m3；r為空間徑向坐標(biāo)，m。

假設(shè)芯塊中心線處的徑向溫度梯度為零，則徑向軸對稱條件和表面溫度邊界條件為：

(2)

式中：rfi和rfo分別為芯塊中心線處半徑和芯塊外表面半徑；Tfo為芯塊外表面溫度。

溫度分布計(jì)算模型的關(guān)鍵之一是材料熱導(dǎo)率的計(jì)算[10]。芯塊及包殼熱導(dǎo)率可通過材料物性分析模型得到。芯塊-包殼間隙傳熱系數(shù)考慮由氣體導(dǎo)熱、輻射傳熱及接觸導(dǎo)熱3種組成，即：

hgap=hsolid+hr+hgas

(3)

式中：hgap為間隙傳熱系數(shù)；hsolid為接觸傳熱系數(shù)；hr為輻射傳熱系數(shù)；hgas為氣體傳熱系數(shù)。

1.3 變形計(jì)算模型

燃料應(yīng)力應(yīng)變計(jì)算是包殼失效分析的基礎(chǔ)，且會(huì)通過芯塊-包殼間隙尺寸直接影響溫度分布計(jì)算。特定溫度、內(nèi)外壓力作用下燃料變形物理過程描述方程包括平衡方程、幾何方程以及本構(gòu)方程[11]。

平衡方程為：

(4)

幾何方程為：

(5)

本構(gòu)方程為：

(6)

(7)

(8)

在外力作用下，材料可能發(fā)生屈服，塑性變形材料應(yīng)變的函數(shù)可表示為σ=f(ε)，如圖2所示。則對于塑性應(yīng)變?chǔ)臥，屈服應(yīng)力可通過下式求得：

圖2 應(yīng)力-應(yīng)變曲線Fig.2 Stress-strain curve

(9)

可寫為：

(10)

該非線性方程可通過牛頓迭代法求解得到：

(11)

1.4 裂變氣體釋放及內(nèi)壓計(jì)算模型

伴隨燃料在堆內(nèi)的裂變反應(yīng)會(huì)產(chǎn)生氣態(tài)和固態(tài)裂變產(chǎn)物。裂變產(chǎn)物會(huì)引起燃料腫脹，加重芯塊與包殼的機(jī)械相互作用(PCMI)。氣態(tài)裂變產(chǎn)物在燃料內(nèi)溶解度較低，會(huì)從燃料芯塊釋放到燃料棒氣空間，增加燃料棒內(nèi)壓，限制燃料棒的壽命；同時(shí)降低棒內(nèi)氣體熱導(dǎo)率，從而使燃料溫度升高。故裂變氣體釋放行為一直都是燃料元件性能分析的重要內(nèi)容之一。

裂變氣體累積產(chǎn)量與發(fā)生的裂變次數(shù)相關(guān)：

(12)

式中：GPT為該區(qū)域裂變氣體累計(jì)產(chǎn)量；Bu為該區(qū)域累積燃耗；VF為該區(qū)域燃料體積；Av為阿伏伽德羅常數(shù)；PR為裂變氣體生成比率，主要考慮的裂變氣體為Kr和Xe。

裂變氣體通常認(rèn)為通過擴(kuò)散進(jìn)行釋放，利用晶內(nèi)擴(kuò)散及晶間行為模擬。假設(shè)晶體為球體，裂變氣體從中心點(diǎn)向外擴(kuò)散，通過求解擴(kuò)散方程可得到晶界位置的裂變氣體濃度，基本擴(kuò)散方程[12]為：

(13)

2 程序開發(fā)

2.1 程序結(jié)構(gòu)

Athena程序系統(tǒng)設(shè)計(jì)如圖3所示。為開展全堆芯反應(yīng)堆物理、熱工水力、燃料性能耦合分析，設(shè)計(jì)接口包括輸入輸出接口、物理耦合接口、熱工耦合接口。輸入輸出接口用于向燃料軟件提供初始輸入及計(jì)算結(jié)果輸出。物理耦合接口與反應(yīng)堆物理計(jì)算軟件耦合，為燃料提供功率分布。熱工耦合接口與熱工水力計(jì)算軟件耦合，為溫度模塊提供邊界溫度分布。程序耦合接口可與CVR 1.0耦合，更加真實(shí)地描述反應(yīng)堆中燃料棒工況，為全堆芯并行化分析提供基礎(chǔ)。

圖3 燃料元件性能分析程序系統(tǒng)Fig.3 Fuel performance analysis code system

2.2 多物理耦合方案

Athena程序多物理耦合計(jì)算流程如圖4所示。通過中子物理模塊獲得燃料功率分布，經(jīng)溫度模塊計(jì)算出燃料棒溫度分布，利用力學(xué)模塊得到燃料應(yīng)力-應(yīng)變更新燃料包殼間隙尺寸，通過間隙傳熱系數(shù)進(jìn)行熱工-力學(xué)耦合；隨后計(jì)算裂變氣體釋放及燃料氣腔內(nèi)壓，更新內(nèi)壓并反饋到溫度和力學(xué)模塊計(jì)算中，當(dāng)氣壓達(dá)到收斂時(shí)完成耦合計(jì)算。

圖4 程序耦合計(jì)算流程圖Fig.4 Coupling flow chart of code

2.3 程序并行化

程序并行方式如圖5所示。利用程序接口定義詳細(xì)燃料功率、冷卻劑邊界條件后燃料棒之間相互獨(dú)立[13]，Athena程序通過消息傳遞接口(MPI)技術(shù)中的單程序多數(shù)據(jù)(SPMD)方法實(shí)現(xiàn)多棒并行分析，使用主進(jìn)程將讀取的輸入文件廣播給所有子進(jìn)程，以便所有子進(jìn)程都能獲取到輸入文件中的參數(shù)信息，并根據(jù)讀取的參數(shù)信息進(jìn)行對應(yīng)燃料棒溫度、應(yīng)力應(yīng)變、內(nèi)壓、包殼腐蝕等行為的數(shù)值計(jì)算，進(jìn)程執(zhí)行完成后，得到所有棒的性能分析結(jié)果，實(shí)現(xiàn)反應(yīng)堆燃料元件性能并行化處理。

圖5 燃料元件多棒并行方式Fig.5 Multi fuel rod parallel mode

3 程序計(jì)算結(jié)果

采用典型商用壓水堆核電站燃料數(shù)據(jù)和同類軟件計(jì)算結(jié)果對程序進(jìn)行初步驗(yàn)證。對標(biāo)同類軟件為法國原子能委員會(huì)(CEA)、法國電力公司(EDF)和法瑪通公司聯(lián)合開發(fā)的輕水堆燃料元件性能分析程序Meteor。

3.1 溫度場分析

燃料棒平均燃耗和功率隨時(shí)間的變化如圖6所示。燃料棒所在的組件在反應(yīng)堆內(nèi)經(jīng)歷兩個(gè)輻照循環(huán)，第1個(gè)循環(huán)長度為394 d，第2個(gè)循環(huán)長度為526 d，兩個(gè)循環(huán)燃料組件在堆芯內(nèi)的位置不同，平均功率略有不同。燃料棒燃耗隨輻照時(shí)間加長而逐漸加深。燃料芯塊峰值溫度隨時(shí)間的變化如圖7所示。反應(yīng)堆啟動(dòng)后，由于芯塊“重定位”作用，芯塊直徑會(huì)瞬間增大，間隙減小，間隙熱導(dǎo)隨之增加，芯塊溫度降低。隨著反應(yīng)堆的運(yùn)行，芯塊和包殼間隙逐漸減小甚至接觸，即芯塊-包殼間隙閉合。間隙閉合后，間隙傳熱系數(shù)相對穩(wěn)定，燃料芯塊中心溫度變化與功率變化相似，壽期末，隨著燃耗加深，燃料芯塊熱導(dǎo)率降低，雖然功率降低但芯塊峰值溫度稍有升高。

圖6 燃料棒平均燃耗及平均線功率隨時(shí)間的變化Fig.6 Fuel rod average burnup and average power vs time

圖7 芯塊峰值溫度隨時(shí)間的變化Fig.7 Fuel pellet peak temperature vs time

3.2 變形計(jì)算

壽期末燃料芯塊及包殼直徑沿軸向分布如圖8所示。燃料包殼受到內(nèi)、外壓作用，由于冷卻劑外壓大于內(nèi)壓，在長期蠕變作用下包殼壽期末直徑小于制造直徑。由于受到燃料元件兩端端塞影響，燃料包殼兩端蠕變變形較小，而分析程序未考慮端塞的支撐作用，同時(shí)兩端芯塊腫脹、熱膨脹等變形較小，造成包殼兩端向內(nèi)蠕變變形較大。對于燃料元件中間區(qū)段，兩程序計(jì)算結(jié)果與測量值符合較好。

圖8 燃料元件壽期末直徑沿軸向分布Fig.8 Axial distribution of fuel cladding diameter at end of life

燃料元件包殼軸向伸長隨時(shí)間的變化如圖9所示。由圖9可看出，在壽期初包殼由于蠕變等作用長度縮短，隨燃耗逐漸增加，輻照生長變形明顯增加，包殼軸向變長[14]。平均直徑、燃料元件包殼伸長數(shù)據(jù)列于表1。

表1 壽期末燃料元件尺寸結(jié)果Table 1 Fuel element dimension result at end of life

圖9 包殼軸向伸長隨時(shí)間的變化Fig.9 Axial elongation of cladding vs time

3.3 裂變氣體釋放

裂變氣體釋放隨時(shí)間的變化如圖10所示，裂變氣體釋放份額隨燃耗的增加逐漸增加，當(dāng)氣態(tài)裂變產(chǎn)物在晶界內(nèi)累積且燃料溫度高時(shí)，裂變氣體釋放速率明顯增加。

圖10 裂變氣體釋放隨時(shí)間的變化Fig.10 Fission gas release vs time

裂變氣體釋放數(shù)據(jù)列于表2。試驗(yàn)單獨(dú)測量了Kr和Xe兩種裂變氣體元素釋放，其中Xe釋放量大于Kr。Athena相較Meteor程序，燃料棒內(nèi)腔體積及內(nèi)壓計(jì)算結(jié)果與測量結(jié)果吻合更好，總體裂變氣體釋放數(shù)據(jù)預(yù)測合理。

表2 裂變氣體釋放數(shù)據(jù)Table 2 Data of fission gas release

4 結(jié)論

本工作開發(fā)了數(shù)值反應(yīng)堆原型系統(tǒng)CVR1.0燃料性能分析程序Athena，能夠?qū)崿F(xiàn)燃料元件熱工-力學(xué)行為高性能模擬，具備多棒并行計(jì)算能力。利用典型商用壓水堆核電站數(shù)據(jù)及同類程序計(jì)算結(jié)果對Athena程序進(jìn)行了初步驗(yàn)證，對比分析了燃料溫度、變形、裂變氣體釋放及內(nèi)壓等計(jì)算結(jié)果，結(jié)果表明Athena程序計(jì)算結(jié)果可靠。

開發(fā)高精細(xì)、強(qiáng)耦合、高效并行運(yùn)行的軟件系統(tǒng)是準(zhǔn)確模擬核反應(yīng)堆燃料性能的基礎(chǔ)[15]。CVR1.0旨在充分利用國產(chǎn)超算的優(yōu)勢解決核反應(yīng)堆復(fù)雜的工程問題，最終將結(jié)合多尺度、多物理、先進(jìn)并行技術(shù)推動(dòng)核反應(yīng)堆分析軟件創(chuàng)新發(fā)展。Athena程序充分利用并行技術(shù)在精細(xì)建模中的優(yōu)勢，為CVR1.0全堆芯多物理耦合提供了重要支持。