包覆燃料顆粒應力的有限元分析

曹 彬，劉 兵，唐春和

（清華大學 核能與新能源技術研究院，北京 100084）

包覆燃料顆粒應力的有限元分析

曹 彬，劉 兵，唐春和

（清華大學 核能與新能源技術研究院，北京 100084）

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

包覆燃料顆粒；應力；有限元法；非球形

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

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

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

1 包覆燃料顆粒SiC層應力計算

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

1.1 計算幾何模型

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

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

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

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

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

1.2 應力計算方法

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

1.3 有限元應力分析

1）有限元分析模型描述

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

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

2）有限元計算參數變化

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

表1 參數變化Table 1 Parametric range

2 結果分析及討論

2.1 球形包覆燃料顆粒應力分析

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

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

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

圖2 球形顆粒最大應力與SiC楊氏模量、SiC層厚度、SiC內徑和內壓的關系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

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

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

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

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

2.2 橢球形包覆燃料顆粒應力分析

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

圖5 橢球顆粒最大應力與球形度的關系Fig.5 Max stress of ellipsoidal particle for variation in eccentricity

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

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

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

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

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

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

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

2.3 凸臺形包覆燃料顆粒應力計算

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

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

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

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

圖11示出球形度1.04、內壓17MPa、R2＝160μm時應力在凸臺形包覆燃料顆粒中的分布。應力在凸臺形包覆燃料顆粒中的分布相較于球形及橢球形包覆燃料顆粒復雜很多?？煽吹?，應力在兩個部分較大，一是過渡圓弧的內側，二是平臺封頭的外側中心處。一般，最大應力出現在平臺的外側中心處，第1段圓弧處的應力與相同內壓下球形SiC層中的應力基本一致。

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

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

3 結論

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

1）有限元法可方便有效地分析包覆燃料顆粒的應力。在后續工作中，可使用有限元法對更多的破損因素進行考慮。有限元法的主要優點是易對二維或三維的因素進行處理，簡化力學推導的過程。

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

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

［1］ 唐春和.高溫氣冷堆燃料元件［M］.北京：化學工業出版社，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］.原子能科學技術，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］ 楊林，劉兵，邵友林，等.高溫氣冷堆橢球形包覆燃料顆粒應力及破損率［J］.清華大學學報：自然科學版，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）；清華大學自主科研計劃資助項目（20121088038）

曹 彬（1990—），男，山東滕州人，碩士研究生，核科學與技術專業