包覆燃料顆粒應(yīng)力的有限元分析

曹 彬，劉 兵，唐春和

（清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084）

包覆燃料顆粒應(yīng)力的有限元分析

曹 彬，劉 兵，唐春和

（清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084）

高溫氣冷堆的燃料元件由包覆燃料顆粒彌散在石墨基體中組成。在反應(yīng)堆運行過程中，輻照及各復(fù)雜的物理化學(xué)反應(yīng)產(chǎn)生的應(yīng)力會使包覆燃料顆粒發(fā)生破損，對包覆燃料顆粒進(jìn)行應(yīng)力分析是評價燃料元件和反應(yīng)堆運行安全性能的主要內(nèi)容之一。本文基于壓力殼模式，主要考慮內(nèi)壓作用下的球形殼層應(yīng)力及包覆燃料顆粒的非球形因素，用有限元法對應(yīng)力進(jìn)行了分析。

包覆燃料顆粒；應(yīng)力；有限元法；非球形

燃料元件是反應(yīng)堆的核心部件，直接影響反應(yīng)堆的經(jīng)濟(jì)性和安全性。燃料元件設(shè)計的目標(biāo)是最有效地生產(chǎn)和導(dǎo)出裂變能，在整個壽期內(nèi)保持包覆材料的完整性，最大限度地約束燃料和放射性裂變產(chǎn)物，保證工作人員和周邊環(huán)境的安全。高溫氣冷堆球形燃料元件由直徑約50mm的燃料區(qū)和厚約5mm的無燃料區(qū)組成，包覆燃料顆粒均勻彌散在燃料區(qū)內(nèi)。TRISO型包覆燃料顆粒在直徑約0.5mm的UO2核芯表面有4層包覆層。從里往外第1層是低密度熱解碳，第2層是高密度各向同性熱解碳，第3層是SiC，第4層也是高密度各向同性熱解碳。整個包覆燃料顆粒的直徑約1mm［1］。

包覆燃料顆粒的破損會造成放射性裂變產(chǎn)物的不可接受性釋放，同時影響反應(yīng)堆的可靠運行。造成包覆燃料顆粒破損的機制有很多，許多國家為適應(yīng)特定的需要研究了自己的包覆燃料顆粒破損模型，這些模型考慮了不同的破損因素。德國的FZJ模型［2］和日本的JAERI模型［3］均主要考慮壓力殼式破損：FZJ模型只考慮SiC層，未考慮熱解碳層帶來的效果，JAERI模型考慮了熱解碳層的收縮和膨脹及每層的單獨失效機制。英國的STRESS3模型［4］著眼于包覆層的應(yīng)力分析，考慮了裂變產(chǎn)物和CO2氣體、核心腫脹、元件制造時核芯燒結(jié)產(chǎn)生的殘留應(yīng)力、熱膨脹系數(shù)的不匹配、輻照帶來的尺寸變化及顆粒和基體之間的相互作用。美國的PARFUME模型［5］考慮了幾乎所有已知的TRISO包覆燃料顆粒破損機制，它借助于有限元軟件對諸如外熱解碳層的破裂和非球形等三維因素進(jìn)行分析。法國的ATLAS模型［6］也考慮了幾乎所有已知的燃料元件堆內(nèi)行為和破損機制，它的一個特點是明確考慮了燃耗大于10%時燃料的結(jié)構(gòu)變化。俄羅斯的GOLT模型［7］正在發(fā)展中，最新的GOLT-v3模型考慮了幾乎所有的輻照、熱機械、化學(xué)以及熱傳遞因素，在正常工況和事故工況下均能充分模擬燃料元件的行為。清華大學(xué)初步建立了包覆燃料顆粒計算模型和球形石墨燃料元件計算模型［8］，主要考慮壓力殼式破損。

包覆燃料顆粒計算模型的建立對包覆燃料顆粒的制造及輻照實驗的設(shè)計起積極指導(dǎo)作用，但目前尚無公認(rèn)的計算模型。模型中考慮的破損因素日漸完善，對非球形因素的討論卻很少。我國的計算模型在不斷發(fā)展，但還未用數(shù)值方法對一些破損因素進(jìn)行考慮。本文用有限元法對包覆燃料顆粒在正常工況下的非球形破損因素進(jìn)行計算。

1 包覆燃料顆粒SiC層應(yīng)力計算

在包覆燃料顆粒的4層包覆層中，SiC層是主要的承壓層。因此，在初步計算中，主要考慮單層SiC在內(nèi)壓作用下的應(yīng)力分布。包覆燃料顆粒在正常制備條件下是理想或近乎理想的球形顆粒，但由于制造誤差，一些非球形的包覆燃料顆粒被壓制在燃料元件中，這些包覆燃料顆粒的應(yīng)力分布與球形包覆燃料顆粒的應(yīng)力分布有很大區(qū)別，對包覆燃料顆粒的應(yīng)力進(jìn)行計算是最重要的研究內(nèi)容［9］。

1.1 計算幾何模型

在實際的包覆燃料顆粒中，非球形的幾何形狀具有多樣性，但有兩種幾何結(jié)構(gòu)很具代表性，分別是橢球形包覆燃料顆粒和凸臺形包覆燃料顆粒，其他包覆燃料顆粒的形狀均可由這兩種顆粒形狀綜合而成。為減少計算量，不使用復(fù)雜的三維計算模型，而采用相對簡單的二維軸對稱的計算模型。在幾何建模中，采用球形度ε對非球形程度進(jìn)行描述。

橢球形包覆顆粒的幾何模型如圖1所示，橢圓內(nèi)徑的長軸長為2a，短軸長為2b，厚度為d。將此圖形繞對稱軸旋轉(zhuǎn)1周即可得到三維的橢球殼。該模型中，球形度為長軸長與短軸長之比，即ε＝a／b。

凸臺形包覆燃料顆粒的模型由兩部分球殼及一部分平板封頭平滑連接而成，幾何模型亦示于圖1。其中，大球殼的內(nèi)徑為R1，對應(yīng)圓心角為ψ1；過渡球殼的內(nèi)徑為R2，對應(yīng)圓心角為ψ2；平板封頭的半徑為R3。三者厚度均為d。將此圖形繞對稱軸旋轉(zhuǎn)1周即可得到三維的凸臺形包覆燃料顆粒模型。該模型中，球形度為：

圖1 橢球形和凸臺形包覆顆粒的幾何模型Fig.1 Ellipsoidal and faceted particle geometry

模型中，R1為已知量，選取一定的ε及R2即可完全確定凸臺的形狀［10］。

1.2 應(yīng)力計算方法

在包覆燃料顆粒應(yīng)力-應(yīng)變-位移關(guān)系的計算中，最常采用的方法是解析解法，該方法將問題簡化為一維球，忽略了非球形因素、溫度梯度和其他一些因素，使計算時間大為縮短［2］。另一種方法是有限元法，它能考慮更多導(dǎo)致包覆燃料顆粒破裂的因素，避免了解析解法中的一些簡化，但需更大的計算量。有限元法在解決球形度問題上有很大優(yōu)勢：1）能清晰得到SiC球殼中每點的應(yīng)力，直觀反映位移以及應(yīng)力的變化趨勢；2）能簡化力學(xué)推導(dǎo)的復(fù)雜過程，以便將更多精力用到結(jié)果分析及更多破損因素的考慮上；3）能分析其他形狀的非球形包覆燃料顆粒的應(yīng)力。因此，本文選取有限元法對包覆燃料顆粒的應(yīng)力進(jìn)行分析，同時便于后續(xù)考慮更多破損因素時的計算。COMSOL Multiphysics［11］是一款通用的有限元分析軟件，本文使用此款軟件對包覆燃料顆粒的應(yīng)力進(jìn)行分析。

1.3 有限元應(yīng)力分析

1）有限元分析模型描述

非球形包覆燃料顆粒在COMSOL中被描述為二維軸對稱結(jié)構(gòu)，故邊界條件也要設(shè)置為對稱的。SiC是一種黏彈性材料，但其性質(zhì)與線彈性材料的非常接近，在計算中，將SiC設(shè)定為線彈性材料。SiC的材料屬性也是計算中必不可少的，其楊氏模量為：

其中：B和T0為表示材料特性的常數(shù)，B＝0.04GPa／K，T0＝962K；E0為0K時的楊氏模量，一般取室溫下的楊氏模量460GPa。SiC的泊松比隨著雜質(zhì)程度、孔隙率、結(jié)晶度的不同有很大的變化。高溫氣冷堆包覆燃料顆粒所用的SiC是純凈的化學(xué)氣相沉積SiC，一般認(rèn)為在燃料溫度范圍內(nèi)其泊松比是常數(shù)0.21［12］。外部邊界條件設(shè)定為自由的，在后續(xù)考慮熱解碳層的輻照收縮效應(yīng)后，外部邊界條件會變得較復(fù)雜。在有限元網(wǎng)格劃分過程中，隨著球形度的不同，元素的個數(shù)也不同。需保證在網(wǎng)格密度增大的情況下得到的結(jié)果變化非常小，這樣得到的應(yīng)力才可靠。

2）有限元計算參數(shù)變化

根據(jù)包覆燃料顆粒制造的實際情況，很多參數(shù)均會圍繞標(biāo)準(zhǔn)值有一微小變化，這些參數(shù)的變化會對包覆燃料顆粒的應(yīng)力分布產(chǎn)生影響。包覆燃料顆粒中，SiC層所受的內(nèi)壓主要來源于CO、CO2（一般忽略不計）和裂變氣體產(chǎn)物，它的大小與核芯尺寸、重金屬燃耗等一系列堆工參數(shù)有關(guān)。在計算中，內(nèi)壓選取反應(yīng)堆正常運行工況下的參考范圍10～25MPa。計算時主要考慮內(nèi)壓、球形度、SiC層半徑、SiC層厚度的變化，具體參數(shù)列于表1。

表1 參數(shù)變化Table 1 Parametric range

2 結(jié)果分析及討論

2.1 球形包覆燃料顆粒應(yīng)力分析

SiC層的楊氏模量隨溫度的變化而變化，通過計算可得包覆燃料顆粒在反應(yīng)堆正常運行溫度下的楊氏模量變化范圍為410～510GPa。楊氏模量對應(yīng)力分布的影響很微弱（圖2a），完全可忽略，所以在后續(xù)的計算中，SiC的楊氏模量皆設(shè)為460GPa。

圖2b示出了SiC層厚度對包覆燃料顆粒中最大應(yīng)力的影響。球形包覆燃料顆粒內(nèi)壓取17MPa，內(nèi)徑為385μm；橢球形包覆燃料顆粒內(nèi)壓為17MPa，球形度為1.1；凸臺形包覆燃料顆粒球內(nèi)壓取17MPa，球形度為1.04，R2為240μm。對于球形和橢球形包覆燃料顆粒，目前的設(shè)計參數(shù)中SiC層厚度為35μm，減小5μm引起的應(yīng)力增加值在20MPa左右，而增加5μm引起的應(yīng)力減小值在15MPa左右，減小SiC層厚度對應(yīng)力的影響更大一些。在SiC層厚度為35μm時，最大應(yīng)力對許用應(yīng)力來說是安全的，增加厚度對提高安全性的意義不大，而減小SiC層厚度提高的應(yīng)力值較大，因此目前的設(shè)計標(biāo)準(zhǔn)是較理想的。在厚度小于標(biāo)準(zhǔn)值時，凸臺形包覆燃料顆粒的應(yīng)力增加很大，而增加厚度對減小SiC層中的應(yīng)力意義不大。

SiC層內(nèi)徑對球形包覆燃料顆粒應(yīng)力的影響示于圖2c，內(nèi)壓為17MPa、SiC厚度為35μm?？煽闯觯SSiC層內(nèi)徑的增大，最大應(yīng)力的增加值有限，內(nèi)徑增加10μm，最大應(yīng)力增加3MPa左右。因此，對于SiC層的參數(shù)設(shè)計，半徑的變化對應(yīng)力的影響不是很大。但在實際生產(chǎn)中，增大SiC層半徑意味著增大核芯體積，從而增大在反應(yīng)堆正常運行時產(chǎn)生的內(nèi)壓，這會引起應(yīng)力的增加。

圖2 球形顆粒最大應(yīng)力與SiC楊氏模量、SiC層厚度、SiC內(nèi)徑和內(nèi)壓的關(guān)系Fig.2 Max stress of spherical particle for variation in Young’s modulus of SiC，thickness of SiC，inner radius of SiC and internal pressure

內(nèi)壓對球形包覆燃料顆粒中應(yīng)力的影響示于圖2d。在其他條件不變的情況下，隨內(nèi)壓的增大，SiC層中的最大應(yīng)力也逐漸增加，內(nèi)壓每增加4MPa，最大應(yīng)力均增加25MPa左右，但變化梯度基本保持不變?？梢?，內(nèi)壓對SiC層中應(yīng)力的影響是最大的，因此，對內(nèi)壓產(chǎn)生的機制進(jìn)行詳細(xì)研究也很有必要。

圖3、4為內(nèi)壓17MPa、SiC層厚度35μm、內(nèi)徑385μm時得到的應(yīng)力分布，SiC層中的應(yīng)力沿徑向由內(nèi)向外均勻減小，最大應(yīng)力與最小應(yīng)力相差約30MPa。

圖3 應(yīng)力在球形顆粒SiC層中的分布Fig.3 Distribution of stress in spherical SiC layer

圖4 球形顆?？v軸處的應(yīng)力Fig.4 Stress of spherical particle on ordinate axis

2.2 橢球形包覆燃料顆粒應(yīng)力分析

圖5示出內(nèi)壓分別為10、17、25MPa時橢球形包覆燃料顆粒中最大應(yīng)力隨球形度的變化。由于球形度的改變，由內(nèi)壓引起的殼層應(yīng)力有較大變化，在殼體的局部產(chǎn)生了高于理想球形包覆燃料顆粒的應(yīng)力，隨球形度的增加，最大應(yīng)力隨之增加。球形度相同時，內(nèi)壓對最大應(yīng)力的影響很大。

圖5 橢球顆粒最大應(yīng)力與球形度的關(guān)系Fig.5 Max stress of ellipsoidal particle for variation in eccentricity

圖6為內(nèi)壓17MPa、球形度1.2、SiC層厚度35μm時得到的橢球形顆粒應(yīng)力分布。與球形SiC層相同，應(yīng)力沿徑向由內(nèi)向外逐漸減小，不過短軸處的應(yīng)力集中程度明顯高于長軸處的，最大應(yīng)力出現(xiàn)在橢球短軸的內(nèi)徑頂點處。圖7為半長軸和半短軸處應(yīng)力沿坐標(biāo)的變化曲線，可見，應(yīng)力的變化較均勻，變化率與球形包覆燃料顆粒的變化率相似。

圖6 橢球形顆粒應(yīng)力在SiC層中的分布Fig.6 Distribution of stress in ellipsoidal SiC layer

圖7 橢球顆粒半長軸和半短軸處的應(yīng)力Fig.7 Stress of ellipsoidal particle on semi-major and semi-minor axes

最大應(yīng)力與平均應(yīng)力之比可定量表示應(yīng)力集中的程度，其中，平均應(yīng)力由應(yīng)力對體積的積分除以總體積得到。橢球顆粒不同球形度的應(yīng)力集中程度示于圖8。由圖可知，橢球形包覆燃料顆粒的應(yīng)力集中特點不很明顯，應(yīng)力分布較均勻，最大應(yīng)力約為平均應(yīng)力的1.1～1.2倍。應(yīng)力集中程度隨球形度的增加逐漸增大。

圖8 橢球顆粒不同球形度的應(yīng)力集中程度Fig.8 Ratio of peak to average stress of ellipsoidal particle for variation in eccentricity

根據(jù)薄殼理論的拉梅解答和無矩理論也可對橢球形包覆燃料顆粒SiC層的應(yīng)力進(jìn)行計算。對比解析解法與有限元法得到的結(jié)果可驗證計算的有效性。表2列出在相同條件下分別使用兩種方法得到的應(yīng)力。從表2可看出，兩種方法得到的結(jié)果相差約2%，基本一致。

表2 解析解法與有限元法結(jié)果的對比Table 2 Comparison of analytical and finite element method calculations

2.3 凸臺形包覆燃料顆粒應(yīng)力計算

由凸臺形包覆燃料顆粒幾何模型建立的分析可知，凸臺形顆粒的形狀可由球形度和過渡球殼的內(nèi)徑R2決定。內(nèi)壓取17MPa時，凸臺形包覆燃料顆粒的最大應(yīng)力隨球形度和過渡球殼內(nèi)徑的變化規(guī)律如圖9所示。球形度一定時，最大應(yīng)力隨R2的增大逐漸減小，球形度越大，最大應(yīng)力減小的趨勢越明顯。R2一定時，最大應(yīng)力隨球形度的增加非常明顯。在實際生產(chǎn)制造過程中，凸臺形包覆燃料顆粒的過渡球殼半徑約為240μm，球形度超過1.04的凸臺形包覆燃料顆粒的數(shù)量也很少。當(dāng)球形度小于1.04、R2為 240μm左右時，SiC層中的最大應(yīng)力基本在許用應(yīng)力范圍內(nèi)。當(dāng)內(nèi)壓取其他值時，最大應(yīng)力隨球形度及R2變化的規(guī)律相同，只是隨內(nèi)壓的增大，凸臺形包覆燃料顆粒的應(yīng)力增大很多。

圖9 凸臺形顆粒最大應(yīng)力隨R2和球形度的變化Fig.9 Max stress of faceted particle for variation in R2and eccentricity

凸臺形包覆燃料顆粒的應(yīng)力集中程度如圖10所示，選取內(nèi)壓為17MPa、R2為40μm的情況進(jìn)行分析，隨球形度的增加，最大應(yīng)力與平均應(yīng)力之比增加較明顯，當(dāng)球形度大于1.1時，最大應(yīng)力均在平均應(yīng)力的5倍以上。說明凸臺形包覆燃料顆粒的應(yīng)力集中特點很明顯，因此，在實際生產(chǎn)過程中，應(yīng)嚴(yán)格控制球形度較大的凸臺形包覆燃料顆粒的出現(xiàn)。

圖10 凸臺形顆粒不同球形度的應(yīng)力集中程度Fig.10 Ratio of peak to average stress of faceted particle for variation in eccentricity

圖11示出球形度1.04、內(nèi)壓17MPa、R2＝160μm時應(yīng)力在凸臺形包覆燃料顆粒中的分布。應(yīng)力在凸臺形包覆燃料顆粒中的分布相較于球形及橢球形包覆燃料顆粒復(fù)雜很多。可看到，應(yīng)力在兩個部分較大，一是過渡圓弧的內(nèi)側(cè)，二是平臺封頭的外側(cè)中心處。一般，最大應(yīng)力出現(xiàn)在平臺的外側(cè)中心處，第1段圓弧處的應(yīng)力與相同內(nèi)壓下球形SiC層中的應(yīng)力基本一致。

圖11 凸臺形顆粒應(yīng)力在SiC層中的分布Fig.11 Distribution of stress in faceted SiC layer

凸臺形包覆燃料顆粒SiC層的應(yīng)力計算無其他詳細(xì)結(jié)果可對比。但得出的凸臺形包覆燃料顆粒的應(yīng)力集中特點明顯，對破損率有較大影響，最大應(yīng)力出現(xiàn)在平臺處與Miller等［10］的結(jié)論基本一致。

3 結(jié)論

基于壓力殼模式，主要考慮了內(nèi)壓作用下的球形殼層應(yīng)力及包覆燃料顆粒的非球形因素，用有限元法對應(yīng)力進(jìn)行了分析。根據(jù)分析和討論，得到以下結(jié)論：

1）有限元法可方便有效地分析包覆燃料顆粒的應(yīng)力。在后續(xù)工作中，可使用有限元法對更多的破損因素進(jìn)行考慮。有限元法的主要優(yōu)點是易對二維或三維的因素進(jìn)行處理，簡化力學(xué)推導(dǎo)的過程。

2）橢球形包覆燃料顆粒和凸臺形包覆燃料顆粒的應(yīng)力分布有很大不同，這兩種非球形顆粒具有代表意義，其他形狀的非球形包覆燃料顆粒均可由這兩種顆粒組合而來。

3）當(dāng)考慮外致密熱解碳層的輻照尺寸變化時，在正常工況下，收縮的熱解碳會給SiC層一壓應(yīng)力，減小SiC層中的應(yīng)力。本文將SiC層作為唯一的承壓層，在這種假設(shè)下，計算得到的應(yīng)力偏大，但對包覆燃料顆粒的設(shè)計制造有一定的指導(dǎo)意義。凸臺形包覆燃料顆粒的應(yīng)力集中程度明顯大于橢球形和球形包覆燃料顆粒的，尤其在球形度較大時，因此，在包覆燃料顆粒的制作過程中應(yīng)盡量減少球形度較大的凸臺形包覆燃料顆粒的出現(xiàn)。

［1］ 唐春和.高溫氣冷堆燃料元件［M］.北京：化學(xué)工業(yè)出版社，2007：33-36.

［2］ LIU Bing，LIANG Tongxiang，TANG Chunhe.A review of TRISO-coated particle nuclear fuel performance models［J］.Rare Metals，2006，25（6）：337-342.

［3］ SAWA K，SHIOZAWA S，MINATO K，et al.Development of coated particle failure model under high burnup irradiation［J］.Journal of Nuclear Science and Technology，1996，33（9）：712-720.

［4］ MARTIN D G.Considerations pertaining to the achievement of high burn-ups in HTR fuel［J］.Nuclear Engineering and Design，2002，213（2）：241-258.

［5］ MILLER G K，PETTI D A，VARACALLE D J，Jr，et al.Statistical approach and benchmarking for modeling of multi-dimensional behavior in TRISO-coated fuel particles［J］.Nuclear Materials，2003（317）：69-82.

［6］ POWERS J J，WIRTH B D.A review of TRISO fuel performance models［J］.Journal of Nuclear Materials，2010，405（1）：74-82.

［7］ GOLUBEV I E，KURBAKOV S D，SCHERNIKOV A.Experimental and computational study of the pyrocarbon and silicon carbide barriers of HTGR fuel particle［J］.Atomic Energy，2008，105（1）：18-31.

［8］ 楊林，劉兵，邵友林，等.HTR-10球形燃料元件模型分析［J］.原子能科學(xué)技術(shù)，2007，41（增刊）：326-330.

YANG Lin，LIU Bing，SHAO Youlin，et al.Program of fuel element in HTR-10［J］.Atomic Energy Science and Technology，2007，41（Suppl.）：326-330（in Chinese）.

［9］ 楊林，劉兵，邵友林，等.高溫氣冷堆橢球形包覆燃料顆粒應(yīng)力及破損率［J］.清華大學(xué)學(xué)報：自然科學(xué)版，2010，50（6）：891-895.

YANG Lin，LIU Bing，SHAO Youlin，et al.Stress and failure probability of ellipsoidcoated particle fuels for an HTGR［J］.Journal of Tsinghua University：Science and Technology，2010，50（6）：891-895（in Chinese）.

［10］MILLER G K，WADSWORTH D C.Treating asphericity in fuel particle pressure vessel modeling［J］.Journal of Nuclear Materials，1994，211（1）：57-69.

［11］馬惠，王剛.Comsol Multiphysics基本操作指南和常見問題解答［M］.北京：人民交通出版社，2009：3-10.

［12］SNEAD L L，NOZAWA T，KATOH Y，et al.Handbook of SiC properties for fuel performance modeling［J］.Journal of Nuclear Materials，2007，371（1-3）：329-377.

Stress Analysis of Coated Particle Fuel Using Finite Element Method

CAO Bin，LIU Bing，TANG Chun-he
（Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing100084，China）

The fuel element of high temperature gas-cooled reactor is composed of coated particle fuel which is dispersed in graphite matrix.In normal operation，the stress due to irradiation and a variety of complex physical and chemical reactions will cause failure of the coated particle fuel.Therefore，the stress analysis of coated particle fuel is important for the safety of fuel element and reactor.The stress was analyzed by the finite element method based on the inner pressure failure mechanism considering asphericity of the particles.

coated particle fuel；stress；finite element method；asphericity

TL352.28

A

1000-6931（2014）02-0318-07

10.7538／yzk.2014.48.02.0318

2012-11-17；

2012-12-27

國家科技重大專項資助（ZX06901）；清華大學(xué)自主科研計劃資助項目（20121088038）

曹 彬（1990—），男，山東滕州人，碩士研究生，核科學(xué)與技術(shù)專業(yè)