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

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

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

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

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

1 數(shù)理模型

1.1 程序介紹

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

1.2 溫度分析模型

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

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

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

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

(1)

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

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

(2)

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

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

hgap=hsolid+hr+hgas

(3)

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

1.3 變形計算模型

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

平衡方程為：

(4)

幾何方程為：

(5)

本構(gòu)方程為：

(6)

(7)

(8)

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

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

(9)

可寫為：

(10)

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

(11)

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

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

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

(12)

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

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

(13)

2 程序開發(fā)

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

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

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

2.2 多物理耦合方案

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

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

2.3 程序并行化

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

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

3 程序計算結(jié)果

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

3.1 溫度場分析

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

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

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

3.2 變形計算

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

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

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

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

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

3.3 裂變氣體釋放

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

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

裂變氣體釋放數(shù)據(jù)列于表2。試驗單獨測量了Kr和Xe兩種裂變氣體元素釋放，其中Xe釋放量大于Kr。Athena相較Meteor程序，燃料棒內(nèi)腔體積及內(nèi)壓計算結(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é)行為高性能模擬，具備多棒并行計算能力。利用典型商用壓水堆核電站數(shù)據(jù)及同類程序計算結(jié)果對Athena程序進行了初步驗證，對比分析了燃料溫度、變形、裂變氣體釋放及內(nèi)壓等計算結(jié)果，結(jié)果表明Athena程序計算結(jié)果可靠。

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