鈉冷快堆燃料元件性能分析程序的開發與驗證

陳啟董,高付海

(中國原子能科學研究院 核工程設計研究所,北京 102413)

從安全角度看,精確預測快堆燃料元件在輻照環境下的演變越來越重要。在燃料元件運行期間,輻照引起的燃料和包殼的熱和力學性能的變化非常復雜,因其行為不但取決于制造過程,而且取決于功率水平、功率歷史、燃耗等運行因素[1]。預測或評估使用壽命期間的燃料行為的一種有效方法是使用燃料性能分析程序進行分析計算。鈉冷快堆燃料元件性能分析程序FIBER從2011年開始開發,用于分析處理正常和預期瞬態條件下的鈉冷快堆燃料元件的行為,程序不僅可為輻照實驗的解釋提供重要信息,還可為安全分析、輻照實驗前的探索性分析、模型驗證等提供重要信息。

FIBER程序代碼由許多理論模型、經驗模型以及控制計算過程的參數組成。然而,燃料行為不能僅通過這些模型的簡單組合來解釋,因為燃料的行為是許多現象相互耦合的結果。因此,必須考慮幾種類型的模型組合來充分再現燃料行為。從這個角度來看,必須使用盡可能多的實際測量案例來進行代碼驗證,以確定合適的模型及參數選擇。

本文首先對FIBER程序的主要模型進行介紹,然后將調研得到的俄羅斯示范快堆BN600燃料元件的實驗數據與FIBER程序計算結果進行比較,并在裂變氣體釋放和芯塊包殼間隙、輻照變形等方面進行討論。

1 燃料元件性能分析的計算模型

1.1 燃料元件總體介紹

FIBER程序是中國原子能科學研究院開發的鈉冷快堆燃料元件性能分析程序,最初版本1.0為分析中國實驗快堆燃料元件(最大燃耗6at%)的行為而開發,為滿足示范快堆的應用(最大燃耗10at%),對程序功能及模型進行了擴展及驗證,目前程序的版本為2.0。程序整體上可分為熱分析模型、力學分析模型。程序的熱學分析主要包括冷卻劑傳熱、燃料和包殼的導熱、燃料-包殼間隙的傳熱、裂變氣體釋放、重結構、裂變產物遷移(MOX)等模型。力學分析主要包括彈性、塑性(包含熱壓)、蠕變、燃料和包殼的接觸應力、芯塊開裂和重定位等模型。程序包含的主要材料特性如下:1) MOX燃料(PuO2含量范圍從0～30%)或UO2燃料;2) 奧氏體不銹鋼(316(Ti)SS、15-15Ti)和鐵素體-馬氏體鋼(HT-9,EP-450)。

燃料元件的計算模型在軸向上分為若干段,每個軸向段將在徑向上分成一些同心環。對于熱和力學分析,都是基于軸對稱的基本假設。

圖1為FIBER程序的分析流程圖。對于每個時間步的分析流程,首先通過線功率、冷卻劑溫度、間隙傳熱系數來進行燃料與包殼溫度的計算。根據溫度計算結果,計算裂變氣體的釋放與裂變產物腫脹,計算芯塊重結構和裂變產物遷移。然后進行彈性、塑性、蠕變、接觸應力的力學分析。熱分析與力學分析以時間步的形式進行推進。通過時間步長的控制確保程序的穩定運行。

圖1 FIBER程序的分析流程圖Fig.1 Analysis flow chart of FIBER code

1.2 熱分析模型

在燃料元件的熱分析中,一維導熱模型用于燃料元件徑向溫度計算,溫度的計算采用有限體積法,冷卻劑溫度是熱分析的邊界條件,在穩定和瞬態條件下,每個軸向節點的冷卻劑溫度可以用功率和冷卻劑流量計算得到,芯塊與包殼的溫度通過導熱計算得到。在這個導熱模型中1個芯塊一般被分成100個節點,包殼在徑向上一般被分成10個節點,以確保有足夠的計算精度。

1) 熱導率模型

對于MOX燃料的熱導率模型,使用以下模型來評價混合氧化物燃料的孔隙率和密度的影響[2]。

λ=K1dK1pK2pK3xK4rλ0

K2p=1-αpP

K3x=1

(1)

式中:λ為熱導率,W/(m·K);OM為氧金屬比;Bu為燃耗,MW·d/kgU;αp為氣孔修正系數;P為氣孔率;T為芯塊溫度,K。

2) 間隙傳熱模型

芯塊與包殼之間的間隙傳熱模型采用了改進的Ross和Stoute模型[3-5]。該模型由3部分組成:氦氣與裂變氣體的氣體導熱;包殼和芯塊表面的固體接觸傳熱;包殼和芯塊表面之間的輻射傳熱。模型可評估在穩態時隨著燃耗而變化的間隙寬度下的間隙傳熱,以及在瞬態時隨著接觸力的變化而變化的固體接觸傳熱。間隙傳熱模型如下:

(2)

式中:hgas、hcon、hr分別為氣體、接觸、輻射換熱系數,W/(cm2·K);λgas為混合氣體的熱導率,W/(cm·K);δ為間隙,cm;C為粗糙度修正系數;R1、R2分別為芯塊、包殼表面粗糙度,cm;g1、g2為氣體躍遷距離,cm;λm為芯塊與包殼等效熱導率,W/(cm·K);pcon為芯塊與包殼接觸壓力,MPa;R為等效粗糙度,cm;H為包殼的硬度,MPa;r1、r2分別為芯塊與包殼尺寸,cm;ε1、ε2為輻射發射率;σ為斯蒂芬玻爾茲曼常數;T1、T2分別為芯塊、包殼溫度,K。

1.3 芯塊重結構模型

因為快堆燃料芯塊徑向功率的分布是均勻的,其溫度高于壓水堆燃料的,在快堆燃料中,會發生由孔隙遷移和晶粒生長引起的重結構現象(中心孔和柱狀晶粒形成)[6]。中心孔的模擬是通過氣孔遷移模型[7]實現的,柱狀晶區的形成是通過晶粒長大模型與氣孔遷移實現的。

P(r,0)=P0

(3)

式中:r為半徑,cm;vp為氣孔遷移速率,cm/s;Δr為網格寬度,cm;Δt為時間步長,s;P0為初始氣孔率。

通過質量守恒即可計算得到中心孔尺寸:

(4)

1) 氣孔遷移的速率

孔隙的遷移速率考慮了溫度與溫度梯度的影響[8-9],具體如下:

(5)

2) 晶粒長大模型

晶粒的長大模型考慮了溫度與裂變氣體的影響[10],有:

(6)

K=5.24×107exp(-2.67×105/RT)

(7)

1.4 裂變氣體釋放模型及輻照腫脹模型

1) 裂變氣體釋放模型

整個裂變氣體釋放模型的計算流程如圖2所示,首先計算晶粒的生長尺寸,根據晶粒生長后的尺寸插值得到新的裂變氣體分布。然后根據裂變率計算晶粒內部氣泡的尺寸、氣泡的密度。通過求解球形晶粒的擴散方程,計算得到擴散進入晶粒邊界的裂變氣體數量。根據溫度與外界環境條件,計算晶粒邊界氣泡可儲存的最大裂變氣體量,裂變氣體的量大于閾值時,裂變氣體將被釋放到氣腔。

圖2 裂變氣體釋放計算流程圖Fig.2 Fission gas release flow chart

裂變氣體擴散方程模型基于理想球形晶粒的方程給出[8,11],裂變產生的氣體原子溶于基體中,由于濃度梯度擴散的作用向晶粒邊界移動。溶于基體的氣體原子與儲存在晶粒內部氣泡的氣體原子的數量取決于捕獲與再溶解的平衡。氣體原子在晶粒中的擴散和被晶粒內部氣泡捕獲由以下方程描述:

(8)

式中:c為裂變氣體晶粒基體中的原子數,cm-3;D為氣體原子的擴散系數,cm2/s;g為氣體原子被晶粒內部氣泡捕獲的速率,s-1;b′為氣體重新溶于晶粒基體的速率,s-1;m為氣泡內單位體積原子數,cm-3;q為單位體積芯塊氣體原子產生率,cm-3·s-1。

晶粒內部的氣泡數量通過WHITE模型[11]進行計算,氣泡的尺寸變化是通過晶內氣泡的氣體狀態方程計算得到:

(9)

(10)

裂變氣體擴散到達晶界,在晶界上形成了棱鏡狀的氣泡[12],棱鏡狀的氣泡所能儲存的裂變氣體的數量是通過氣體狀態方程來確定的。當晶界氣泡儲存的量超過閾值時,氣體將形成隧道,氣體立即從晶粒釋放到氣腔。

(11)

式中:pgas為氣體壓力,Pa;R為氣泡半徑,cm;n為裂變氣體數量;f(θ)為晶粒邊界氣泡形狀函數。

2) 輻照腫脹

裂變氣體導致的輻照腫脹由晶粒內部氣泡與晶界氣泡兩部分組成。

(12)

晶粒內部氣泡體積腫脹為:

(13)

晶界氣泡體積腫脹為:

(14)

式中:K為氣泡的數量;Vgrain為晶粒體積,cm3;f(θ)為棱鏡狀氣泡的形狀函數,右邊分母中的2表示1個氣泡是由2個晶粒共享的。

3) 氣腔壓力

計算中假設裂變氣體是理想氣體,并且棒中的壓力是均勻的,計算方法如下:

(15)

1.5 力學分析模型

力學分析模型中,應力/應變分析采用有限元方法,用四自由度的四邊形單元進行,每個節點都有一個自由度,芯塊徑向的網格劃分數可以是10～100。假設所有節點的軸向位移相同。使用虛功原理,時間tn+1的平衡條件[13]表示如下:

(16)

作用于燃料元件的芯塊和包殼的外力如圖3所示,在軸向上,包殼受到彈簧產生的軸向力,包殼外表面受到一回路冷卻劑產生的壓力,內表面除了受到裂變氣體產生的壓力,還受到芯塊對包殼的接觸作用力。

圖3 作用于芯塊和包殼的力Fig.3 Forces acting on fuel and cladding

力學分析計算的流程如圖4所示。首先計算與應力無關的變形,例如致密化、輻照腫脹、熱膨脹等,然后計算單元的彈性模量等材料屬性,然后組裝剛度矩陣,通過剛度矩陣計算應力-應變,應變收斂后通過迭代確保接觸應力的計算是收斂的。假如出現接觸到無接觸,無接觸到接觸狀態的改變的現象,應采用時間步截斷的方式,在芯塊-包殼間隙為0、接觸應力為0的時間點截斷,所有與時間相關的物理量根據時間進行插值,重新設計時間步長進行下一步的計算。

圖4 力學計算流程圖Fig.4 Mechanical flow chart

1.6 時間步控制

為保證程序計算過程的穩定,采用了如下的時間步控制。

1) 在1個時間步長內,線功率的變化在10 W/cm以內[4]。

2) 在1個時間步長內,燃耗的變化在0.05at%以內[4]。

3) 在裂變氣體釋放的模塊中,時間增量Δt是由擴散方程模型求解方法決定的[4,11,13]。

Δt≤0.05(ΔR)2/D

(17)

式中:ΔR為網格的寬度,cm;D為裂變氣體的擴散系數,cm2/s。

4) 在有限元彈塑性計算中,為保證蠕變計算的穩定,應使蠕變應變增量不超過彈性應變[13]。

(18)

5) 為保證物質遷移的穩定性,應確保網格內物質的含量不會出現負值,即:

Δt≤nVj/(Jj,l-Jj,r)

(19)

式中:J為通過截面的物質的量,mol/s;n為單位體積物質的量,mol/cm3;V為體積,cm3。

6) 在非穩態的溫度計算過程中,時間步長的限制如下[14-15]。

(20)

式中,a為熱擴散系數,cm2/s。

7) 芯塊-包殼接觸應力計算中的時間步截斷[14-15]。

對于時間步執行直到t→t+Δt的計算,芯塊與包殼的接觸狀態從非接觸到接觸變化或從接觸到非接觸式變化時,會對時間步截斷,減小時間步長Δt,確定間隙寬度變為零的時間或接觸壓力變為零的時間t′,所有時間相關的物理量都會隨時間線性插值。然后,改變在時間t′的接觸狀態,并重新計算從t′到t+Δt的物理量。

2 快堆燃料元件的算例驗證

為保證鈉冷快堆燃料元件性能分析程序FIBER選用的模型及計算假設的合理性,應采用實驗數據驗證程序計算的準確性。本文調研得到了俄羅斯示范快堆BN600的UO2燃料元件和MOX燃料元件的尺寸和運行參數,具體列于表1與圖5。BN600-標準燃料元件為示范快堆BN600的驅動燃料元件,該燃料元件的包殼材料為1515Ti系列奧氏體不銹鋼,芯塊為17%富集度的UO2,穩態運行3個周期,累計滿功率559 d。該燃料元件主要進行了裂變氣體釋放量、氣體體積、輻照腫脹等關鍵參數的測量[16-20]。BN600-MOX燃料元件為示范快堆BN600的試驗燃料元件,該燃料元件的包殼材料為1515Ti系列奧氏體不銹鋼,芯塊為MOX燃料,包殼直徑略小于標準燃料元件,燃料元件穩態運行3個周期,累計滿功率559 d。該燃料元件主要進行了柱狀晶區、包殼尺寸、間隙等關鍵參數的測量[16-20]。

表1 BN600燃料元件參數[16-19]Table 1 BN600 fuel element parameter[16-19]

圖5 相對功率分布Fig.5 Relative power distribution

采用鈉冷快堆燃料元件性能分析程序FIBER對表1中的燃料元件建模計算,燃料元件沿徑向分環,沿軸向分段進行計算。

2.1 裂變氣體釋放的驗證

圖6為BN600-標準燃料元件最大燃耗4at%～10at%裂變氣體的測量結果[17,20],對30根燃料元件進行穿刺,得到了標準狀態下的氣體體積(包含初始充入氦氣),圖中橫坐標為標準狀態下的氣體體積除以燃料元件中芯塊的質量,縱坐標為燃料元件最大燃耗。30根燃料元件均在BN600反應堆第一流量區穩態運行,組件流量相同,運行過程中燃料組件并不移動位置,燃料元件的運行時間不同,所達到的最大燃耗不一致。從圖6可看出,FIBER程序計算結果與實驗結果的趨勢一致,可很好地反映氣腔中裂變氣體的變化。

圖6 裂變氣體釋放量計算結果與實驗結果對比Fig.6 Comparison of fission gas release calculation results with experimental data

圖7為7根BN600-標準燃料元件燃耗裂變氣體壓力的測量結果[17,20],氣體的壓力是在20 ℃下測量得到。從圖7可看出,FIBER程序計算結果與實驗結果的趨勢一致,最大相對誤差為26%。

圖7 氣腔壓力計算結果與實驗結果對比Fig.7 Comparison of plenum pressure calculation results with experimental data

2.2 芯塊輻照變形

文獻[17,20]測量得到的UO2芯塊的輻照變形,具體數據如圖8所示,當燃耗超過4at%,燃料元件的腫脹隨燃耗線性增加。在燃耗為9at%～10at%時,體積變形為9%～10%,在這種情況下,UO2芯塊的徑向輻照變形率為每1%燃耗為0.33at%。需要說明的是,原始數據為示范快堆運行的燃料元件的測量結果,未記錄初始燃料芯塊的尺寸,變形率按照表1的名義尺寸進行統計,燃料芯塊的公差為-0.15～0 mm,因此圖中出現了徑向變形率小于0的情況。采用程序對上下公差芯塊進行計算,從圖8可看出,程序計算結果與實驗結果的趨勢一致。

圖8 UO2芯塊變形計算結果與實驗數據的對比Fig.8 Comparison of UO2 pellet strain results with experimental data

大部分的MOX燃料的徑向變形率[17,20]與UO2芯塊的徑向變形率相同,即每1%燃耗為0.33at%,但部分芯塊變形率達到了每1%燃耗為0.6at%～0.7at%。這部分燃料的徑向變化率有一部分裂變產物遷移的貢獻,在輻照后MOX燃料與包殼的間隙檢測到了Mo、Cs等裂變產物形成的氧化物(簡稱JOG=Joint Oxyde-Gaine),如圖9所示。同樣需要說明的是,原始數據為示范快堆運行后燃料元件的測量結果,未記錄初始燃料芯塊的尺寸,變形率按照表1的名義尺寸進行統計。燃料芯塊的公差為-0.15～0 mm,采用程序對上下公差芯塊進行計算,并考慮JOG的貢獻,從圖9可看出,程序計算結果與實驗結果的趨勢一致。

圖9 MOX芯塊變形計算結果與實驗結果的對比Fig.9 Comparison of MOX pellet strain results with experimental data

2.3 包殼尺寸

圖10為MOX試驗件的測量結果[17,20],對兩根燃料元件7個位置進行了測量,橫坐標為相對高度(距離燃料段底端高度/燃料段總高度,下同),縱坐標為包殼管外徑測量結果。包殼的尺寸從原始的6.6 mm最大增加到6.9 mm。需要說明的是,由于周圍棒束及外套管的約束作用,輻照腫脹后包殼管的截面不再是圓形,截面呈現橢圓形,如圖11所示,同一燃料元件位置的數據包含了橢圓形長軸與短軸的測量結果。

圖10 包殼變形計算結果與實驗數據的對比Fig.10 Comparison of cladding strain results with experimental data

圖11 輻照前后包殼形狀Fig.11 Shape of cladding before and after irradiation

從圖10可看出,程序計算結果與實驗結果的趨勢一致。實驗數據與程序計算結果差異較大的點出現在相對位置0.19處,在實驗中該位置觀察到了芯塊較大的變形(每1%燃耗為0.7at%)。該位置包殼的變形可能是芯塊較大的變形引起的接觸應力導致的。程序計算結果與實驗數據的偏差可能是程序只考慮了徑向的物質遷移,未考慮軸向裂變產物遷移,芯塊變形率計算偏小導致的,同時程序假設JOG并不會對芯塊與包殼的接觸有貢獻,只對傳熱有貢獻。

根據輻照數據及程序結算結果推測,活性區中平面附近芯塊的外表面溫度1 050 ℃,間隙溫度755 ℃,高溫下JOG為液態。但隨著時間推移,活性中平面的JOG沿軸向遷移至燃料元件的底端,該位置的溫度相對較低,芯塊的外表面溫度622 ℃,間隙溫度521 ℃,JOG在燃料元件的底端凝固。銫會與燃料發生化學反應[21],引入銫與燃料的化學反應腫脹,該物質會對芯塊與包殼的接觸應力有很大的貢獻。目前程序對JOG的模型不夠完善,后續考慮基于輻照現象對該部分模型進行修改。

2.4 間隙尺寸對比

圖12為MOX試驗件的測量結果[17,20]。由于芯塊開裂以及包殼輻照變形后形狀的不規則,同一位置包含多個測量結果。從圖12可看出,程序計算結果與實驗結果的趨勢一致。相對位置0.93處出現較大的偏差,可能是由于該位置的芯塊初始芯塊尺寸偏小。

圖12 間隙計算結果與實驗結果對比Fig.12 Comparison of gap calculation results with experimental data

2.5 柱狀晶區尺寸

圖13為MOX試驗件的柱狀晶粒范圍的測量結果[17,20]。從圖13可看出,程序計算結果與實驗結果的趨勢一致。

圖13 柱狀晶粒范圍的計算結果與實驗結果對比Fig.13 Comparison of columnar region calculation results with experimental data

3 結論

本文對鈉冷快堆燃料元件性能分析程序FIBER的熱分析、芯塊重結構、裂變氣體釋放、力學、時間步控制模型進行了介紹,并通過調研獲取了俄羅斯BN600反應堆兩種類型燃料元件的輻照后檢驗結果,通過實驗數據驗證了FIBER的計算結果。結果表明,程序可模擬最高燃耗11.8at%、輻照損傷78 dpa燃料元件輻照變形、柱狀晶區、裂變氣體釋放的現象。