楊 爍，呂俊男，李 群,*

(1.西安交通大學 航天航空學院 機械結構強度與振動國家重點實驗室，陜西 西安 710049； 2.中國核動力研究設計院 反應堆燃料及材料重點實驗室，四川 成都 610213)

彌散型燃料元件[1-2]是將易裂變?nèi)剂舷嘁孕☆w粒的形式均勻彌散在相對惰性的非裂變材料基體中的一種混合物核燃料。由于燃料顆粒被基體材料分隔開，使裂變產(chǎn)物損傷局限在燃料顆粒及其周圍基體內(nèi)，不受裂變產(chǎn)物損傷的基體能形成連續(xù)網(wǎng)絡，這樣有利于包容裂變產(chǎn)物、抑制燃料腫脹，且燃料元件芯部溫度低，導熱性好，從而允許燃料達到更高的燃耗和安全可靠性，因此被廣泛用于研究試驗堆和動力堆[3]。

相對于普通的顆粒增強復合材料，彌散燃料的特殊性在于中子輻照條件下，核燃料顆粒內(nèi)的易裂變物質，如鈾會發(fā)生裂變和嬗變，此過程中產(chǎn)生的裂變氣體不斷聚集形成氣泡，并導致燃料顆粒輻照腫脹[2]。高燃耗狀態(tài)下，氣體腫脹致使燃料顆粒逐漸演變成多孔結構，顆粒內(nèi)氣孔的尺寸和位置的非均勻分布與孔間的應力干涉效應極易導致氣孔開裂貫通。事故工況下，燃料芯體內(nèi)部顆粒破裂的比例很高，大量微裂紋貫穿、連通擴展至基體，最終宏觀表現(xiàn)為包殼表面起泡，可認為燃料顆粒的開裂失效是彌散燃料失效的根源[4-5]。

為研究彌散燃料的失效行為，需建立相應的理論模型，發(fā)展有效的數(shù)值模擬方法，學者們在這方面開展了許多有意義的工作。Rest等[6-8]研究了燃料內(nèi)裂變氣體腫脹模型的計算以及U-Mo合金中裂變氣泡尺寸分布規(guī)律與高溫成核機制。萬遠富等[9]利用Eshelby等效夾雜理論，在考慮溫度場和輻照腫脹的條件下，分析了顆粒的形狀以及輻照腫脹對應力狀態(tài)的影響。龍沖生等[10]在裂變氣孔尺寸均勻分布的條件下利用彈性力學建立了彌散燃料顆粒開裂模型，系統(tǒng)計算預測了燃料顆粒的開裂規(guī)律。Jeong等[11-12]建立了彌散燃料的局部應力失效分析模型，并用有限元模型模擬了燃料基體和顆粒的熱-力學行為。Ding等[13-15]采用代表性體元并假設燃料顆粒均勻分布，建立了細觀有限元模型，分析了細觀熱-力學行為和輻照蠕變行為對彌散型板狀燃料元件可靠性的影響。趙毅等[16-18]基于數(shù)值模擬對高燃耗結構的燃料顆粒內(nèi)部氣孔力學行為進行研究，建立了雙氣泡和均勻分布的多氣泡模型，討論了燃料顆粒內(nèi)部氣孔的尺寸、分布等對其應力分布的影響，得到了燃料顆粒裂紋起源是在表層氣泡內(nèi)壁的結論。目前的研究都是基于顆粒或氣孔均勻周期性分布情況下展開的，這與實際不符，開展考慮燃料顆粒內(nèi)部非均勻性的數(shù)值模擬，能更好地反映顆粒的強度特性。

本文基于隨機序列吸附法(RSA)思想實現(xiàn)含彌散非均布的裂變氣孔的高燃耗陶瓷燃料顆粒的幾何建模，并結合ABAQUS有限元軟件二次開發(fā)實現(xiàn)燃料顆粒靜力學分析的參數(shù)化建模。分析關鍵參數(shù)氣孔尺寸、溫度、及顆粒所受基體材料約束壓應力對燃料顆粒內(nèi)部最大拉應力的影響規(guī)律，并對燃料顆粒內(nèi)部危險區(qū)的分布進行討論。

1 燃料顆粒參數(shù)化建模

彌散燃料芯體中的每一個燃料顆粒都可單獨看作1個微型燃料元件，在輻照條件下燃料顆粒內(nèi)部彌散分布了大量的氣孔。在高溫和高燃耗條件下，氣孔內(nèi)部裂變氣體的壓力會達到上百MPa[19]，氣孔的尺寸和位置分布均會改變?nèi)剂项w粒內(nèi)部應力場的分布，進而影響開裂行為?？紤]這些參數(shù)的影響必須參數(shù)化建模，通過參考輻照后高燃耗燃料顆粒微觀形貌的電鏡照片可知，燃料顆粒及其內(nèi)部彌散分布的小氣孔近似為球形[20]。因此，本文對燃料顆粒進行抽象與簡化(圖1)：1) 燃料顆粒為球形，顆粒由致密的燃料相和離散的裂變氣孔兩相組成；2) 燃料顆粒外表面受到均布的約束壓應力pf，源于基體材料對燃料相腫脹變形的約束作用，彌散氣孔內(nèi)表面受到均布的氣體壓應力pg作用。

圖1 高燃耗燃料顆粒的抽象與簡化Fig.1 Abstraction and simplification of high burnup fuel particle

文獻[21]利用對數(shù)正態(tài)分布對燃料顆粒中的氣孔半徑數(shù)據(jù)進行過擬合，但這種方法對不同燃耗的適用范圍有要求。本文采用適應性更強的Weibull分布來描述氣孔尺寸及氣孔間距的分布規(guī)律。Weibull分布的概率密度函數(shù)如下：

(1)

式中：x為隨機變量，本文中代表氣孔半徑或氣孔間距；λ為尺度參數(shù)，反映變量的平均特性；k為形狀參數(shù)，反映變量的非均勻性，k越大則離散性越小。當k取不同值時，可得到正偏、負偏和對稱的概率密度函數(shù)。

參數(shù)化建模流程如圖2所示。首先，基于RSA思想，運用Python編程手段，在燃料顆粒內(nèi)部投放小氣孔，并記錄氣孔的坐標和半徑，此過程保證氣孔相互不重疊，并控制氣孔半徑(r)和相鄰氣孔間距(d)都服從Weibull分布，分別用Weibull-r和Weibull-d表示，以顆粒半徑160 μm、平均氣孔半徑2 μm的模型為例，所繪效果圖如圖3所示。然后，使用Python腳本語言對ABAQUS有限元進行二次開發(fā)，建立顆粒幾何模型，并實現(xiàn)程序化加載、劃分網(wǎng)格和靜力學計算。

圖2 參數(shù)化建模流程Fig.2 Parametric modeling process

圖3 氣孔尺寸和位置非均勻分布的 高燃耗燃料顆粒模擬效果圖Fig.3 Simulation effect picture of high burnup fuel particle with heterogeneous distribution of pores

2 數(shù)值模擬

2.1 燃料顆粒的力學參數(shù)

UO2陶瓷核燃料具有耐高溫、耐輻照、腫脹小、包容裂變產(chǎn)物性能好等優(yōu)點，能達到高燃耗的狀態(tài)，是目前應用最廣泛的陶瓷核燃料。本文針對高燃耗狀態(tài)下的UO2燃料顆粒開展數(shù)值模擬研究，定性分析關鍵參數(shù)對燃料顆粒開裂的影響規(guī)律。作為陶瓷型核燃料，UO2燃料顆粒從輻照初始狀態(tài)到高燃耗斷裂失效，致密的燃料相都處于彈性階段，不考慮其塑性變形行為。目前對于輻照后陶瓷顆粒的力學性能研究較少，本文參考文獻[10]中陶瓷的彈性模量E與燃耗和溫度T的關系式：

E=162 000+63 000/

(1+35BU)-20(T+273)

(2)

其中：BU為燃耗深度，F(xiàn)IMA；T為溫度，℃。

燃料顆粒中裂變氣孔體積分數(shù)ρg與燃耗的關系為：

ρg=0.021+0.66BU

(3)

有限元模型中氣孔的最大尺寸與最小尺寸相差倍數(shù)不大，所以假設顆粒內(nèi)部氣孔的內(nèi)壓大小相同，氣孔內(nèi)壓pg采用超高壓下實際氣體的狀態(tài)方程計算：

(4)

(5)

其中：ng為裂變氣體總量；Vg為氣體總體積；R=8.31 J/(mol·K)為氣體常數(shù)；a=5.57×10-5m3、b=2.39×10-5m3為實際氣體狀態(tài)方程參數(shù)；Df=10.96×106/(1+ρg)為燃料顆粒的密度，g/m3；β為裂變氣體(Xe+Kr)的裂變產(chǎn)額；Mf=269 g/mol為燃料相的摩爾質量。

此外，認為圖1所示裂變氣孔內(nèi)壓作用的致密燃料相球殼內(nèi)壁的拉應力大于UO2燃料相的斷裂強度時，燃料顆粒開裂，即選用第一強度理論作為UO2燃料顆粒的開裂判據(jù)。對應有限元模擬研究中，計算得到燃料顆粒內(nèi)部氣孔相互干涉的第一主應力場，其中最大應力點處即為整個燃料顆粒最大拉應力位置。

2.2 有限元模型

本文選取的高燃耗燃料顆粒尺寸Rf定為160 μm，彌散氣孔尺寸Rg也在實驗觀測范圍內(nèi)。欲建立高溫、高燃耗UO2燃料顆粒的全尺寸仿真計算模型，還需對各氣孔內(nèi)壁和顆粒外表面施加均勻壓力(圖4)，載荷數(shù)值由程序代碼計算并控制施加，運用有限元求解得到應力場。此外，由于本研究有限元模型中考慮了氣孔尺寸和分布位置的非均勻性，致使燃料顆粒內(nèi)部的結構非常復雜，為保證計算結果的精確性，在有限元網(wǎng)格劃分過程中，采用四節(jié)點平面應力縮減積分單元CPS4R，對所有的氣孔孔壁和燃料顆粒外表面進行加密分割，完成網(wǎng)格劃分后，每個模型的單元數(shù)都在20萬以上。

圖4 高燃耗UO2燃料顆粒有限元模型的 加載方式Fig.4 Loading mode of high burnup UO2 fuel particle

3 結果分析與討論

3.1 氣孔尺寸對燃料顆粒最大拉應力的影響

為方便定性分析，參數(shù)化建模過程中保證顆粒內(nèi)部氣孔間距滿足一致的非均勻分布形式，而氣孔尺寸相同，選取的氣孔半徑Rg分別為0.8、1.2、1.6、2.0、2.4 μm，利用有限元法計算不同模型的最大拉應力?；诿商乜_方法[22]，在給定氣孔率和氣孔間距滿足Weibull分布形狀參數(shù)的情況下，對每種氣孔尺寸的顆粒模型建立多個隨機模型，取計算結果的平均值作為有效數(shù)據(jù)。其中，仿真模型的材料參數(shù)由式(2)～(5)計算得到，具體如下：顆粒半徑，160 μm；溫度，600 ℃；氣孔率，20%；基體壓應力，0.1 MPa；彈性模量，156 .4 GPa；泊松比，0.316。

在給定氣孔率ρg、溫度T、燃耗BU和基體壓應力pf的情況下，燃料顆粒內(nèi)部的最大拉應力隨氣孔半徑的變化規(guī)律如圖5所示。由圖5可知，氣孔尺寸增大，顆粒的最大拉應力隨之增大，即當氣孔半徑由0.8 μm增大到2.4 μm時，最大拉應力從222.7 MPa增大到484.5 MPa。此過程是非線性變化的，氣孔尺寸越大，最大拉應力增加的幅度越大。如氣孔尺寸從0.8 μm變化到1.2 μm時，最大拉應力增大34 MPa，氣孔尺寸從2.0 μm變化到2.4 μm時，最大拉應力增加122.1 MPa?；谝陨辖Y果，可認為太小的氣孔對燃料顆粒開裂失效的影響不大，大尺寸氣孔對燃料顆粒開裂失效的影響更顯著。

圖5 氣孔半徑對最大拉應力的影響Fig.5 Influence of pore radius on maximum tensile stress

有限元計算的最大拉應力云圖示于圖6。由圖6可看出，不同氣孔尺寸的燃料顆粒最大拉應力出現(xiàn)在靠近顆粒表層，而從高燃耗燃料顆粒的電鏡照片[5]可知，燃料顆粒也確實多從表層附近開始起裂，本研究結果與輻照實驗觀察結果具有一致性。

圖6 最大拉應力出現(xiàn)位置Fig.6 Position of maximum tensile stress

3.2 基體約束壓應力對燃料顆粒最大拉應力的影響

燃料元件在堆內(nèi)服役時受到的環(huán)境靜水壓力通過基體的傳遞最終作用在燃料顆粒上，此外顆粒輻照腫脹也會受到基體的束縛，所以可將基體材料對燃料顆粒的約束壓應力pf看作燃料顆粒開裂行為的阻力，它與顆粒的強度密切相關[23]。選取氣孔半徑Rg=2.4 μm的燃料顆粒模型控制變量分析最大拉應力的變化趨勢，取pf為1、10、20、30、40、50 MPa進行有限元計算，其余材料和模型參數(shù)不變，研究pf對顆粒強度的影響，結果示于圖7。分析圖7可知，在給定氣孔率ρg、燃耗BU、溫度T、氣孔半徑尺寸Rg的情況下，基體對燃料顆粒的約束壓應力pf越大則燃料內(nèi)部最大拉應力越??；pf由1 MPa增大到50 MPa時，最大拉應力近似線性減小了128.5 MPa。可見基體約束壓應力pf對UO2燃料顆粒的開裂行為影響很大，不可忽略。

圖7 基體約束壓應力對最大拉應力的影響Fig.7 Influence of constrained compressive stress of matrix on maximum tensile stress

3.3 溫度對燃料顆粒最大拉應力的影響

堆內(nèi)運行工況下，燃料顆粒處于高溫狀態(tài)，且溫度是變化的，通常與燃料顆粒在燃料元件中的分布位置和運行工況有關。因此，有必要分析溫度變化對燃料顆粒強度的影響。對于氣孔尺寸和間距均滿足Weibull分布且氣孔平均半徑為2.4 μm的燃料顆粒，在給定ρg為20%、pf為0.1 MPa時，溫度增大，其內(nèi)部最大拉應力的變化規(guī)律如圖8所示。圖8表明，氣孔尺寸Weibull分布參數(shù)λ和k不同時，燃料顆粒內(nèi)部最大拉應力均隨溫度升高而增大，且均近似呈線性關系。從圖8還可知，溫度是燃料顆粒失效的敏感參數(shù)，溫度從300 ℃上升到700 ℃時，最大拉應力增加超過200 MPa。如氣孔尺寸Weibull分布形狀參數(shù)k=7時，最大拉應力隨溫度升高而增大了241.3 MPa。形狀參數(shù)k反映氣孔尺寸分布的離散性，從結果分析可知k越小，相同溫度條件下，最大拉應力越大。原因是形狀參數(shù)k越小，氣孔尺寸分布越離散，這會增大大尺寸氣孔的比例。

圖8 溫度對最大拉應力的影響Fig.8 Influence of temperature on maximum tensile stress

3.4 燃料顆粒內(nèi)高應力區(qū)的分布

在輻照條件下，UO2陶瓷燃料相本身的斷裂強度σf會隨微結構變化[24]。本文定義陶瓷燃料顆粒內(nèi)部超過燃料相斷裂強度的區(qū)域為高應力區(qū)域，即開裂危險區(qū)，裂紋先從這些區(qū)域產(chǎn)生。在給定氣孔率、燃耗深度、溫度、基體約束壓應力、氣孔尺寸和間距都非均勻分布的條件下，改變陶瓷燃料相的斷裂強度，考察高應力區(qū)的分布，結果如圖9所示。分析可知，當氣孔內(nèi)壓一定，陶瓷燃料相的斷裂強度從275 MPa減小到175 MPa時，高應力(紅色)區(qū)域的面積增大；從圖9也可看出高應力區(qū)的演化過程，首先開裂危險區(qū)先從陶瓷燃料顆粒的表層多處出現(xiàn)，然后在顆粒內(nèi)部某些地方也出現(xiàn)開裂危險區(qū)，之后這些區(qū)域相互連通，甚至可將顆粒貫穿。此外，尺寸較大氣孔邊緣的高應力區(qū)面積更大，也更易與其周圍的氣孔連通。綜上所述，在輻照條件下，由于氣孔尺寸和位置非均勻分布，顆粒內(nèi)部應力干涉效應復雜，致使顆粒的危險區(qū)從多處出現(xiàn)，即顆粒是從多處開裂破壞的。

圖9 燃料相斷裂強度不同時高應力區(qū)域的分布Fig.9 Distribution of high stress regions with different fracture strengths of fuel phase

3.5 氣孔分布方式對燃料顆粒最大拉應力的影響

高燃耗結構的氣孔尺寸和位置本身是隨機分布的，因此這種模型孔間產(chǎn)生干涉效應是復雜的，相較于傳統(tǒng)的氣孔均勻分布燃料顆粒模型，二者內(nèi)部的應力場存在差異。本文參照文獻[17]的氣孔分布方式建立氣孔均勻分布的燃料顆粒模型，通過有限元計算得到應力云圖，如圖10所示。對與上述模型具有相同氣孔率、內(nèi)壓載荷、約束條件的氣孔非均分布的顆粒模型進行計算。其中非均勻分布的顆粒模型在給定氣孔分布參數(shù)的情況下隨機生成7次，得到7組計算數(shù)據(jù)，將其與均勻分布的模型進行對比，結果如圖11所示。

圖10 氣孔位置均勻分布的燃料顆粒模型應力云圖Fig.10 Stress contour of fuel particle model with uniform arrangement of pore position

圖11表明，非均勻分布的最大拉應力較均勻分布模型的大，更易開裂，且非均勻分布的燃料顆粒內(nèi)部的高應力區(qū)是多處存在的，更能反映高燃耗燃料顆粒多處開裂的物理現(xiàn)象，而均勻分布的顆粒模型沒有此現(xiàn)象。

圖11 氣孔均勻分布與非均勻分布模型結果對比Fig.11 Comparison of model results of uniform and non-uniform arrangement of pores

4 結論

本文結合Python語言和有限元二次開發(fā)，實現(xiàn)了氣孔尺寸和位置非均勻分布的高燃耗陶瓷燃料顆粒參數(shù)化建模，基于有限元技術對陶瓷燃料顆粒的強度和開裂影響因素進行了分析研究，得到如下主要結論。

1) 在給定氣孔率、燃耗深度、溫度、基體約束壓應力，且僅考慮氣孔位置非均勻隨機分布的條件下，燃料顆粒內(nèi)最大拉應力隨氣孔尺寸的增大而增大，且大尺寸氣孔對最大拉應力的影響更大；在給定氣孔率、燃耗深度、溫度、氣孔尺寸，且僅考慮氣孔位置非均勻隨機分布的條件下，燃料顆粒內(nèi)最大拉應力隨基體約束壓應力的增大近似線性減小。

2) 本研究有限元模型中考慮了溫度這一重要開裂影響參數(shù)，同時模擬計算了氣孔尺寸和非均勻隨機分布下的應力場，發(fā)現(xiàn)燃料顆粒內(nèi)最大拉應力受溫度影響顯著，且隨溫度升高而近似線性增大。

3) 在氣孔尺寸和間距都非均勻分布的條件下，改變陶瓷燃料相的斷裂強度能得到高應力區(qū)的演化過程。顆粒內(nèi)部應力干涉效應復雜，致使多處出現(xiàn)顆粒開裂的危險區(qū)域，即顆粒是從多處開裂破壞的，而氣孔均勻分布的燃料顆粒模型沒有這種現(xiàn)象。

高燃耗燃料顆粒在堆內(nèi)運行環(huán)境下的開裂行為是多因素耦合作用的結果，本研究基于有限元理論結合Python語言，為系統(tǒng)分析這種復雜失效行為提供了技術支持，也為進一步研究整個彌散燃料元件的輻照性能奠定了基礎。