張海青 盧林遠(yuǎn) 王 鵬 嚴(yán) 超 林 俊

（中國科學(xué)院上海應(yīng)用物理研究所 上海201800）

在核反應(yīng)堆發(fā)展歷程中，二氧化鈾（UO2）核燃料由于具有熔點高（3 138.15 K）、非裂變組合元素氧的熱中子俘獲截面低（＜0.000 2 b）、在化學(xué)和結(jié)構(gòu)方面具有較高的輻照穩(wěn)定性，且制造費用較低等優(yōu)勢［1］，是輕水反應(yīng)堆核電站廣泛應(yīng)用的核燃料。UO2作為核燃料的顯著缺點是導(dǎo)熱性能較差（7.7 W·m-1·K-1，300 K）［2］，僅為金屬鈾的十幾分之一，這就降低了燃料芯塊到冷卻劑間的能量輸出效率。這意味著它的工作溫度相當(dāng)高，約2 273.15 K。2011年日本福島核電站事故揭露了UO2核燃料低熱導(dǎo)率引發(fā)的核安全隱患［3］。在福島核事故后，為提高核反應(yīng)堆中燃料元件抵抗嚴(yán)重事故的能力，美國能源部發(fā)起了“研制和開發(fā)事故容錯核燃料（Accident Tolerant Fuel，ATF）計劃”［3］，這種燃料與標(biāo)準(zhǔn)的輕水反應(yīng)堆UO2燃料鋯合金包殼相比，在保持或提高反應(yīng)堆正常運行工況下燃料性能的同時，當(dāng)堆芯主動冷卻系統(tǒng)功能喪失后，維持冷卻狀態(tài)的時間更長，安全裕度更高。

中國在響應(yīng)國家安全高效發(fā)展核電的政策下，也啟動了相關(guān)的研發(fā)計劃，旨在通過研發(fā)事故容錯UO2核燃料來推動新一代核燃料元件的開發(fā)與應(yīng)用。此類燃料是通過在UO2燃料基體中添加高熱導(dǎo)率的物質(zhì)，如氧化鈹（BeO）［4］、碳納米管［5］或碳化硅（SiC）纖維［6］等，來提高其熱導(dǎo)率，降低其工作溫度。此類燃料體系的優(yōu)勢在于對工業(yè)制造體系的改動小、制造成本低。

目前，國內(nèi)外關(guān)于事故容錯熱導(dǎo)率增強(qiáng)型UO2核燃料的研究關(guān)注點主要在于制備工藝的優(yōu)化和導(dǎo)熱性能的改善。而微觀結(jié)構(gòu)的調(diào)控是改善導(dǎo)熱性能的主要手段之一。Pan 等［7］研究了短纖維填充橡膠復(fù)合材料中相鄰兩根短纖維之間的相互作用對應(yīng)力分布產(chǎn)生的影響；金泉等［8］主要對顆粒和纖維混合填充橡膠材料進(jìn)行了研究，得出了平行于加載方向的增強(qiáng)體受了很大的應(yīng)力，隨著纖維和加載方向之間夾角的增大應(yīng)力減小；文獻(xiàn)［9-10］研究表明，粒子間距是決定空間填料分布的一個重要參數(shù)，隨著粒子間距離的增大，導(dǎo)熱系數(shù)呈指數(shù)衰減；文獻(xiàn)［11］研究了橡膠基體的填料粒子對復(fù)合材料熱導(dǎo)率的影響。

實際上，復(fù)合材料的有效導(dǎo)熱系數(shù)由許多因素決定，如顆粒的類型和形狀、顆粒的排列方式、體積大小以及填料與基體之間的接觸熱阻。一般來說，用單純的實驗方法很難區(qū)分這些因素的影響。因此，模擬計算對于指導(dǎo)更好的熱應(yīng)用復(fù)合材料的設(shè)計和制造至關(guān)重要。本文將利用商用有限元軟件Ansys Workbench，采用三維有限元分析的方法，通過模擬填充粒子的復(fù)合材料體積單元在特定溫度邊界條件下的溫度、熱流分布情況，分析在UO2燃料中添加BeO 纖維后，填料的幾何形態(tài)、空間位置取向及填料空間堆積方式、顆粒大小等對熱導(dǎo)率的影響。為事故容錯燃料制備工藝的優(yōu)化和導(dǎo)熱性能的改善提供理論依據(jù)。

1 模型建立

1.1 熱傳導(dǎo)方程

對于各向同性材料，不考慮熱源的穩(wěn)態(tài)熱傳導(dǎo)方程為[12]：

式中：k 為熱導(dǎo)率，又叫導(dǎo)熱系數(shù)，熱傳導(dǎo)系數(shù)，W·m-1·K-1；T為溫度，K。熱傳導(dǎo)問題的邊界條件有三種：

式中：h為熱交換系數(shù)；Te為環(huán)境溫度。式（2）表示邊界的溫度值已知；式（3）表示邊界的外法向熱流密度值已知；式（4）表示熱交換邊界條件。

1.2 熱導(dǎo)率

未經(jīng)輻照、致密的燃料顆粒熱導(dǎo)率k0為［13］：

綜合考慮燃耗、孔隙度和裂變產(chǎn)物的影響，輻照下的顆粒熱導(dǎo)率為：

式中：FD和FP分別為溶解和沉淀的裂變產(chǎn)物的影響因子；FM為孔隙度影響因子；FR為輻照效應(yīng)影響因子；B為燃耗（已裂變原子數(shù)與總可裂變原子數(shù)的比值），%；P為孔隙度，%；s為孔的形狀因子。

燃料顆粒和基體的熱導(dǎo)率［14］為：

1.3 有限元模型

UO2顆粒在基體中的分布，可以近似為均勻分布，有限元模型可以將均勻分布設(shè)置成周期性分布的模型，即只需選取一個代表性體積單元來研究。根據(jù)均勻化理論，選取等效體積單元（Representative Volume Element，RVE）［6］，進(jìn)而得到復(fù)合材料的均勻化參數(shù)。本文采用的是Ansys Workbench 穩(wěn)態(tài)熱分析模塊，對于顆粒填充體系進(jìn)行導(dǎo)熱模擬研究。為提高其模擬的精確度選擇3D模型最佳。本文基體材料均為UO2，填料粒子均為BeO 的復(fù)合材料，其中UO2和BeO 的熱導(dǎo)率見文獻(xiàn)［15］單胞上下面各施加不同的溫度載荷，四周絕熱。

因為模型中含有曲率的曲面，所以網(wǎng)格劃分的總體尺寸控制采用曲率尺寸函數(shù)Curvature，并且經(jīng)過檢測網(wǎng)格質(zhì)量在Element Quality規(guī)則下處在好水平范圍內(nèi)。

用Ansys 對上述穩(wěn)態(tài)傳熱問題進(jìn)行求解后，可讀出熱流密度和熱傳導(dǎo)方向上模型兩端的溫度，從而計算出模型的熱導(dǎo)率；

根據(jù)傅里葉穩(wěn)態(tài)熱傳導(dǎo)規(guī)律，通過模型的熱流密度q為［11］：

式中：keff為模型的等效熱導(dǎo)率；Δz為單元體邊長；ΔT為熱傳導(dǎo)方向上的溫差。因而模型的熱導(dǎo)率keff為：

因為多數(shù)燃料在制作過程中填料顆粒是近似均勻的隨機(jī)排列在基體中的，為了驗證周期排列的適用性，結(jié)合MATLAB軟件，編程生成隨機(jī)小球，將包含隨機(jī)小球的幾何圖形導(dǎo)入Ansys Workbench中，施加和周期排列情況相同的邊界條件，運用同樣的熱導(dǎo)率算法，計算了隨機(jī)排列小球的熱導(dǎo)率，取多次隨機(jī)模型的平均值，并和周期排列結(jié)果做了比較。

2 結(jié)果與討論

2.1 填料的幾何形狀對事故容錯燃料熱導(dǎo)率的影響

為了研究BeO 填料幾何形狀對事故容錯燃料導(dǎo)熱性能的影響，本文分別選取了球體、正方體、不同長徑比的圓柱體，以及與圓柱體長度相等、兩端面為正方形的長方體作為填料分別分布在尺寸為100 μm×100 μm×100 μm的基體中，并且根據(jù)填料在UO2基體中的體積分?jǐn)?shù)為5%分別算出球狀填料的半徑、正方體的邊長及不同長徑比圓柱體的直徑和長度，再利用Ansys Workbench軟件模擬穩(wěn)態(tài)傳熱過程。初始條件：頂面和底面施加相應(yīng)的溫度，單元體四周絕熱。參考中國科學(xué)院上海應(yīng)用物理研究所熔鹽堆運行溫度大概是924.15 K，計算了該溫度時材料的熱導(dǎo)率。

表1給出了BeO顆粒占UO2基體中的體積分?jǐn)?shù)為5%時，事故容錯燃料在常溫300 K時熱導(dǎo)率隨填料形狀的變化，其中D:L為圓柱體直徑與長度的比。從表1可以看出，填料的形狀對熱導(dǎo)率有顯著影響。其中圓柱體填料長徑比不同時，熱導(dǎo)率也不同。圓柱體1直徑為100 μm，及和基體單元邊長相同，圓柱體2～8，圓柱體直徑逐漸減小，長度逐漸增加，到圓柱體9，直徑最小，長度最大，為100 μm。在這個過程中，熱導(dǎo)率先下降后增加，與文獻(xiàn)［11］結(jié)論稍有不同，是文獻(xiàn)［8］的補(bǔ)充。圓柱體1 和圓柱體9 兩種極限情況，熱導(dǎo)率都比球體、正方體及和與圓柱體長度相等的長方體要大，所以事故容錯燃料的BeO填料顆粒，建議選用直徑和長度比大于等于3 或者接近1/3.96的圓柱體。

2.2 圓柱體BeO填料纖維的空間取向和熱導(dǎo)率的關(guān)系

為了驗證圓柱體BeO 填料纖維的空間取向和熱導(dǎo)率的關(guān)系，分別計算了§2.1 中圓柱體1、3、4、5、6、9在溫度為300 K、BeO填料顆粒旋轉(zhuǎn)0°、30°、45°、60°、80°、90°（以圓柱軸的方向垂直于溫度加載方向為旋轉(zhuǎn)起始點，旋轉(zhuǎn)到平行于溫度加載方向時，為旋轉(zhuǎn)90°）時的熱導(dǎo)率，結(jié)果見圖1。

由圖1 可以看出，當(dāng)圓柱體直徑與長度的比為1.5:1 時，旋轉(zhuǎn)角度幾乎不影響（如圖1 中圓柱4），但對于直徑與長度的比值大于1.5:1的圓柱體（圓柱體1、3），起始點熱導(dǎo)率最高，隨著旋轉(zhuǎn)角度的變大，熱導(dǎo)率逐漸減小，當(dāng)軸向與溫度加載方向（z軸）一致時，熱導(dǎo)率最低。而對于直徑與長度的比值小于1.5:1的圓柱體（圓柱體5、6、9），情況恰好相反：隨著旋轉(zhuǎn)角度的增大，熱導(dǎo)率逐漸變大，當(dāng)軸向與溫度加載方向（z軸）一致時，熱導(dǎo)率最大，導(dǎo)熱性能最好。

綜合表1和圖1可以看出，當(dāng)圓柱體軸長和基體單元邊長相等時，及D:L=1:3.96，并且軸的方向和溫度加載方向一致時（圖2（d）），熱導(dǎo)率最大。

目前工業(yè)生產(chǎn)的纖維填料在基體中往往是隨機(jī)取向的，由模擬值發(fā)現(xiàn)纖維取向?qū)w系的導(dǎo)熱性能影響顯著。因此為了追求良好的導(dǎo)熱性能，制造過程中就應(yīng)該使纖維的直徑小于長度，并且軸向取向與熱流傳導(dǎo)方向一致或者夾角盡量小。

表1 300 K時ATF的熱導(dǎo)率隨填料顆粒形狀的變化Table 1 The thermal conductivity of the ATF with different shape of packing at 300 K

2.3 球體BeO填料顆粒的排列方式和體積大小對熱導(dǎo)率的影響

2.3.1 BeO填料顆粒的排列方式對熱導(dǎo)率的影響

我們分別采用了MATLAB、COMSOL 軟件和Ansys Workbench 相結(jié)合的隨機(jī)排列方式和單獨使用Ansys Workbench的規(guī)則排列方式，隨機(jī)排列方式的邊長為50 μm，包覆顆粒體積占比為4.8%，共有30個球，半徑為3.63 μm，見圖2（a）；規(guī)則排列的方法是：1）確定一個基本的柵元；2）用該柵元按照一定的方式填充給定的體積。常用到的基本幾何柵元可以是簡單立方體柵元（Simple Cubic lattice，SC）、體立方體柵元（Body Centred Cubic lattice，BCC）、面立方體柵元（Face Centred Cubic lattice，F(xiàn)CC）（表2）。本文選取的基本柵元中球的半徑和和隨機(jī)排列模型的相同，為3.63 μm，包覆顆粒體積占比為4.8%。三種填充方式的有限元模型見圖2（b～d）。圖3給出不同堆積方式下的熱導(dǎo)率值隨溫度的變化。

由圖3可見，填充體積份額相同時，在常溫300～1 100 K 的范圍內(nèi)，無論是隨機(jī)排列，還是三種周期排列，等效熱導(dǎo)率的值都非常接近。因此，排列方式對等效熱導(dǎo)率的影響可以忽略。

圖1 熱導(dǎo)率隨圓柱體BeO填料顆粒旋轉(zhuǎn)角度的變化（右邊為左邊方框處的放大圖）Fig.1 Changes of thermal conductivity with the rotation angle of cylindrical BeO packing particles(the picture on the right is an enlarged one of which on the left)

表2 基本柵元的特性和柵元內(nèi)球數(shù)Table 2 Characteristics of basic lattice and number of balls

圖2 三種填充方式的有限元模型 (a)隨機(jī)排列，(b)SC排列，(c)BCC排列，(d)FCC排列Fig.2 Three different filling formats of finite element model

圖3 300～1 100 K下填料堆積方式對事故容錯燃料熱導(dǎo)率的影響Fig.3 Influence of packing stacking mode on thermal conductivity of ATF(for 300～1 100 K)

2.3.2 BeO填料顆粒的大小對熱導(dǎo)率的影響

為了研究BeO 填料顆粒的體積大小對熱導(dǎo)率的影響，必須考慮接觸熱阻和界面熱阻的影響。

接觸熱阻是由于接觸表面的粗糙度，使有效接觸面積小于表面重疊面積。對于尺寸為幾微米的顆粒，接觸顆粒-顆粒表面的表面粗糙度無關(guān)緊要，對于顆粒-基體接觸，我們可以假設(shè)聚合物基體與顆粒完全接觸。因此，接觸熱阻很小，可以忽略不計。另一方面，界面熱阻是由于粒子和基體的性質(zhì)不匹配造成的。粒子和基體由不同的材料制成，具有不同的密度、德拜溫度、頻率和聲速，在它們之間的邊界處為聲子流動提供阻力。由于我們在基體中有數(shù)千個粒子，這種阻力對聲子的流動非常重要，從而降低了復(fù)合材料的有效導(dǎo)熱系數(shù)。

如果粒子比較大，界面熱阻的影響可以忽略。有時實際上需要在基體中加入更小的粒子，因為當(dāng)顆粒越小，界面越致密，界面熱阻越大。因此，對于固定的界面熱阻，有一個臨界直徑，用d0表示。當(dāng)顆粒直徑大于d0時，顆粒能提高復(fù)合材料的導(dǎo)熱系數(shù)，反之則會降低復(fù)合材料的導(dǎo)熱系數(shù)。所以，需要對顆粒填充復(fù)合材料界面熱阻的影響進(jìn)行分析。

界面熱阻RIF，又叫Kapitza 熱阻，計算方程如下［16］：式中：ρ 為材料密度；?為聲子速度；下標(biāo)p 和m 分別代表顆粒與基體；Cm為基體等壓熱容。

本文分別計算了低體積分?jǐn)?shù)（BeO 填料顆粒體積分?jǐn)?shù)分別為5%的SC、FCC 模型）和高體積分?jǐn)?shù)（BeO 填料顆粒體積分?jǐn)?shù)為52%的SC 模型和74%FCC 模型）的球形BeO 填料顆粒復(fù)合材料，結(jié)果如圖4所示，有效熱導(dǎo)率keff首先隨著R而顯著增加，然后達(dá)到飽和，達(dá)到飽和的點隨BeO填料顆粒體積分?jǐn)?shù)而增大，即對于顆粒體積分?jǐn)?shù)為5%的復(fù)合材料來說，當(dāng)顆粒直徑為100 μm時，達(dá)到飽和，再增大直徑，熱導(dǎo)率將保持不變。對于BeO填料顆粒體積分?jǐn)?shù)為52%的模型，達(dá)到飽和的點為300 μm，對于BeO 填料顆粒體積為74%的模型，達(dá)到飽和的點為1 000 μm。圖4中黑色直線為基體UO2的熱導(dǎo)率，與其他4 條線的交點分別為51 nm、51 nm、50.3 nm、50.3 nm，即BeO 填料顆粒體積分?jǐn)?shù)分別為5%的SC、FCC 模型，d0值均為51 nm；體積分?jǐn)?shù)為52%的SC模型和74%FCC模型，d0值均為50.3 nm，兩組數(shù)值的偏差為1.4%，由此可知BeO填料顆粒體積分?jǐn)?shù)對d0值的影響可以忽略不計。當(dāng)顆粒直徑小于臨界點d0值時，復(fù)合材料的keff低于基體熱導(dǎo)率，在d0到飽和點的范圍內(nèi)，隨著球半徑的增加，keff顯著增加，當(dāng)球半徑大于飽和點后，keff不再增加。

圖4 復(fù)合材料熱導(dǎo)率隨球形填充顆粒半徑的變化Fig.4 The effective thermal conductivity of composites as a function of the radius of spherical filling particles

2.3.3 與理論值的比較

關(guān)于顆粒復(fù)合材料的等效熱導(dǎo)率的公式，前人做了大量的工作，Verma［17］、Maxwell［18］、Babanov［19］、Brailsford［20］、Singh［21］等給出的結(jié)果分別如下，其中?為顆粒的體積百分比含量：

通過式（15）～（19）計算不同顆粒體積含量和熱導(dǎo)率情況下的等效熱導(dǎo)率的理論值，并和有限元的計算值進(jìn)行比較。圖5 為顆粒體積含量不同時，等效熱導(dǎo)率的理論值和計算值的結(jié)果對比。從圖5可以看出，當(dāng)顆粒體積含量較小時，各理論值和計算值的差值相對較小，隨著體積份額的增大，5個理論公式值之間的差距越來越大，其中Babanov、Singh 與計算值的偏差都在20%以上，Brailsford 在?＞20%以后，偏差也大于20%。Verma、Maxwell 在?≤45%時和計算值比較接近，但是大于45%時這兩個公式也不適用了。結(jié)果表明：高填充體積分?jǐn)?shù)BeO的有效導(dǎo)熱系數(shù)對填充量的變化非常敏感。即使是少量的額外填料，熱導(dǎo)率的提高也是顯著的。為了得到更適合顆粒體積含量較大時的公式，對本工作的有限元法（Finite Element Method，F(xiàn)EM）計算結(jié)果進(jìn)行了擬合，得到了UO2基、BeO為填料顆粒的復(fù)合材料的等效熱導(dǎo)率公式如下：

3 結(jié)語

本文在借鑒了國內(nèi)外對事故容錯燃料的設(shè)計和導(dǎo)熱性能研究的基礎(chǔ)上，針對不同的結(jié)構(gòu)設(shè)計方案進(jìn)行了分析比較。具體得出以下幾個方面的結(jié)論：

1）事故容錯燃料的BeO 填料顆粒，建議選用直徑和長度的比大于等于3或者接近1/3.96的圓柱體。

圖5 半徑相同時等效熱導(dǎo)率與體積分?jǐn)?shù)的關(guān)系Fig.5 Relationship between equivalent thermal conductivity and volume fraction with the same radius

2）當(dāng)圓柱體BeO 填料顆粒D:L=1.5 時，方向?qū)釋?dǎo)率幾乎無影響；當(dāng)D:L 大于1.5 時，軸的方向與溫度加載方向夾角越大，熱導(dǎo)率越高，當(dāng)垂直于溫度加載方向時，熱導(dǎo)率最大；當(dāng)D:L小于1.5時，軸的方向與溫度加載方向夾角越小，熱導(dǎo)率越高，當(dāng)平行于溫度加載方向時，熱導(dǎo)率最大。

3）有效熱導(dǎo)率keff首先隨著顆粒直徑R而顯著增加，然后達(dá)到飽和。當(dāng)顆粒直徑小于臨界點d0值時，復(fù)合材料的熱導(dǎo)率低于基體熱導(dǎo)率；在d0到飽和點的范圍內(nèi)，隨著球半徑的增加，熱導(dǎo)率顯著增加；當(dāng)?shù)竭_(dá)飽和點后，隨著BeO 填料顆粒半徑的增加，復(fù)合材料熱導(dǎo)率將保持不變。

4）Verma、Maxwell給出的經(jīng)驗公式在?≤45%時和計算值比較接近，但是大于45%時這兩個公式也不適用了。本文給出了UO2基、BeO 為填料顆粒的復(fù)合材料在? 大于45%時依然適用的等效熱導(dǎo)率公式。