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

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

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

核燃料元件是反應堆系統中的核心部件，是反應堆運行過程中熱量及裂變反應的主要來源，故燃料元件的堆內性能對反應堆經濟性和安全性有重要影響。正常運行條件下，燃料元件設計必須保證其在高溫、高壓、強輻照環境下的完整性。利用計算模擬的方法可以用較低的成本預測燃料元件堆內行為，達到評價元件堆內輻照服役性能、指導燃料元件設計的目的[1]。因此，燃料元件性能分析程序對燃料研發[2]、安全評價具有重要意義[3]。

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

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

1 數理模型

1.1 程序介紹

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

1.2 溫度分析模型

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

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

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

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

(1)

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

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

(2)

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

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

hgap=hsolid+hr+hgas

(3)

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

1.3 變形計算模型

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

平衡方程為：

(4)

幾何方程為：

(5)

本構方程為：

(6)

(7)

(8)

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

圖2 應力-應變曲線Fig.2 Stress-strain curve

(9)

可寫為：

(10)

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

(11)

1.4 裂變氣體釋放及內壓計算模型

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

裂變氣體累積產量與發生的裂變次數相關：

(12)

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

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

(13)

2 程序開發

2.1 程序結構

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

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

2.2 多物理耦合方案

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

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

2.3 程序并行化

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

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

3 程序計算結果

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

3.1 溫度場分析

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

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

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

3.2 變形計算

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

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

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

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

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

3.3 裂變氣體釋放

裂變氣體釋放隨時間的變化如圖10所示，裂變氣體釋放份額隨燃耗的增加逐漸增加，當氣態裂變產物在晶界內累積且燃料溫度高時，裂變氣體釋放速率明顯增加。

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

裂變氣體釋放數據列于表2。試驗單獨測量了Kr和Xe兩種裂變氣體元素釋放，其中Xe釋放量大于Kr。Athena相較Meteor程序，燃料棒內腔體積及內壓計算結果與測量結果吻合更好，總體裂變氣體釋放數據預測合理。

表2 裂變氣體釋放數據Table 2 Data of fission gas release

4 結論

本工作開發了數值反應堆原型系統CVR1.0燃料性能分析程序Athena，能夠實現燃料元件熱工-力學行為高性能模擬，具備多棒并行計算能力。利用典型商用壓水堆核電站數據及同類程序計算結果對Athena程序進行了初步驗證，對比分析了燃料溫度、變形、裂變氣體釋放及內壓等計算結果，結果表明Athena程序計算結果可靠。

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