張明， 厲井鋼， 劉虓瀚， 盧勇， 朱亞楠

(中廣核研究院有限公司反應(yīng)堆工程軟件研究所， 深圳 518000)

燃料棒包殼是防止放射性物質(zhì)外泄的第一道安全屏障，其堆內(nèi)性能直接影響核電廠運行的安全性與經(jīng)濟性。因此，燃料棒包殼的變形機理及其數(shù)值分析方法一直是反應(yīng)堆核燃料研究的重點。

燃料棒包殼一般由鋯合金[1]制成，在壽命初期由于UO2芯塊的密實作用其體積收縮，導(dǎo)致包殼在軸向產(chǎn)生未支撐段，這是包殼蠕變坍塌的先決條件。在反應(yīng)堆內(nèi)高溫、高壓和輻照等耦合作用下，未支撐段包殼會逐漸向內(nèi)蠕變，最終坍塌失效，蠕變是包殼坍塌的根本原因。包殼坍塌會造成局部熱點影響冷卻劑流道，或者使得臨近燃料棒或?qū)蚬軓澢?，影響反?yīng)堆運行安全。中國壓水堆燃料組件設(shè)計標準[2]規(guī)定：包殼要滿足自立條件，且不能發(fā)生蠕變坍塌。由于包殼的服役工況復(fù)雜，難以直接觀測其變形狀態(tài)，因此需要對包殼的變形和應(yīng)力進行建模計算，常見分析方法包括理論解析建模和開發(fā)專用軟件求數(shù)值解。

包殼管的理論建模可追溯到1959年，Hoff等[3]研究了無限長薄圓柱殼在均勻外壓作用下的穩(wěn)態(tài)蠕變行為，但未考慮包殼管的彈塑性變形和瞬態(tài)蠕變。Bargmann[4]通過考慮彈性變形擴展了Hoff的方法，并指出包殼的彈性變形會極大降低包殼壽命。Nishiguchi[5]假設(shè)應(yīng)力沿著壁厚均勻分布，并給出無限長高溫氣冷堆加熱管的蠕變坍塌臨界時間。以上研究只限于模態(tài)二失穩(wěn)，并未考慮失穩(wěn)后的后屈曲行為。隨后，嚴愛民[6]提出了高溫和外壓作用下有限長圓柱殼的多模態(tài)蠕變失穩(wěn)解法；謝旭華[7]開展了無限長夾心殼的后屈曲研究；Xue[8]研究了無限長圓柱薄殼的局部失穩(wěn)現(xiàn)象，并給出了其后屈曲解；盧勇等[9]基于商業(yè)有限元軟件ABAQUS研究了三維燃料棒包殼管的變形和坍塌行為。

雖然理論建?？傻玫桨鼩す苋渥兪Х€(wěn)和后屈曲的解析解，但理論建模只適用于邊界規(guī)則、幾何模型簡單的問題，對于堆內(nèi)變形行為復(fù)雜的包殼管，通常開發(fā)專用數(shù)值軟件計算包殼的堆內(nèi)狀態(tài)，并判斷其是否符合安全標準。常見的包殼變形分析軟件包括BUCKLE[10]、COLAPX[11]、COVE-1[12]、CEPAN[13]和FROCO[14]等，以上軟件采用無限長管假設(shè)(即平面應(yīng)變假設(shè))和截面上的對稱假設(shè)，在蠕變坍塌變形計算和臨界壽命評估方面存在以下不足。

(1)均基于平面應(yīng)變假設(shè)，是一種二維的力學(xué)模型，其幾何模型為無限長管，這在變形分析中會產(chǎn)生較大誤差，雖然COVE-1[12]、CEPAN[13]和FROCO[14]考慮了對有限長管的修正，但計算誤差仍然較大。

(2)COLAPX[11]、COVE-1[12]、CEPAN[13]和FROCO[14]在截面內(nèi)采用對稱假設(shè)，只分析1/4模型，故只能分析圓管變?yōu)闄E圓管的模態(tài)二失穩(wěn)，而實際包殼未支撐段為有限長度，一般會發(fā)生多模態(tài)失穩(wěn)[8]，以上軟件無法準確模擬。

(3)在數(shù)值計算方法方面，BUCKLE[10]、COVE-1[12]和CEPAN[13]采用有限差分法，計算精度依賴于空間和時間步長，大時間步長的計算結(jié)果精度受限；FROCO[14]采用矩陣變換法將計算域離散為許多小結(jié)塊，未考慮結(jié)塊之間的相互作用，難以精確得到包殼的變形和應(yīng)力狀態(tài)。綜上所述，現(xiàn)有計算軟件的基本假設(shè)和計算方法不盡合理，難以精確計算包殼的應(yīng)力和變形，也難以準確預(yù)測包殼的堆內(nèi)狀態(tài)。盡管商業(yè)軟件在包殼管蠕變變形計算中可以被采用[9]，但其仍不如以上專業(yè)軟件應(yīng)用廣泛和計算方便。

針對現(xiàn)存軟件的不足，亟須研究燃料棒包殼的三維變形機理，即建立熱-輻照-力等多物理場耦合作用下的包殼模型，并開發(fā)出高精度的專用計算軟件?，F(xiàn)基于連續(xù)介質(zhì)力學(xué)力學(xué)理論建立包殼管的控制方程，引入相應(yīng)的蠕變本構(gòu)關(guān)系，并開發(fā)基于有限元的包殼管三維變形數(shù)值計算方法。通過與商業(yè)軟件ABAQUS結(jié)果的對比來驗證本軟件的正確性，并計算真實燃料棒包殼的堆內(nèi)變形。

1 方法

1.1 控制方程

直接建立包殼管的三維幾何模型，無需要求包殼管滿足平面應(yīng)力/應(yīng)變假設(shè)。笛卡爾坐標系下一般彈性體的力學(xué)控制方程如下。

(1)平衡方程

(1)

式(1)中：σ為應(yīng)力；x為坐標，i,j=1、2、3分別表示x、y、z方向。式(1)展開為3個方程。

(2)幾何方程

(2)

式(2)中：ε為應(yīng)變；x為坐標；u為位移。式(2)展開為6個方程。

1.2 彈性本構(gòu)關(guān)系

對于彈性體，要想計算出其位移，必須給出彈性本構(gòu)關(guān)系，這里假設(shè)包殼管在彈性階段為線彈性材料，其應(yīng)力應(yīng)變關(guān)系符合廣義胡克定律，即

σij=2Gεij+λεkkδij

(3)

式(3)中：G=E/[2(1-ν)]為剪切模量；λ=Eν/[(1+ν)(1-2ν)]為體積模量，E為彈性模量，ν為泊松比；εkk為體積應(yīng)變。式(3)包含6個方程。

一般地，彈性力學(xué)方程式(1)～式(3)共包含15個方程，需要求15個未知數(shù)，分別為應(yīng)力6個、應(yīng)變6個和位移3個，因此以上3組方程是封閉的方程系統(tǒng)。引入包殼管的蠕變本構(gòu)方程，增加方程的同時也增加了蠕變相關(guān)的未知量，故方程系統(tǒng)還是封閉、適定的。

1.3 蠕變本構(gòu)關(guān)系

(4)

式(4)中：σeff為等效應(yīng)力；T為溫度；φ為中子注量率，其他參數(shù)的取值見表1。參考塑性力學(xué)理論的Prandtl-Reuss流動法則，將等效蠕變應(yīng)變流動到各個方向，即

表1 蠕變本構(gòu)模型中的材料參數(shù)[15]

(5)

式(5)中：Sij為偏應(yīng)力。

1.4 有限元計算方法

以上方程不僅個數(shù)多，而且涉及蠕變引起的材料非線性，無法直接求解，故對以上控制方程進行有限元離散。將有限長的包殼管離散為六面體，采用八節(jié)點線性等參元建模。為了兼顧效率和精度，經(jīng)網(wǎng)格收斂性測試后包殼管徑向、環(huán)向和軸向的網(wǎng)格數(shù)分別為3、40和7個，總網(wǎng)格數(shù)為840，結(jié)構(gòu)化網(wǎng)格劃分如圖1所示。

圖1 包殼管三維網(wǎng)格

八節(jié)點等參元的形函數(shù)為

(6)

式(6)中：ξ、η、ζ分別為母單元上的3個方向的坐標值。

單元內(nèi)任一點的坐標和位移可以用8個節(jié)點的相應(yīng)量插值表示為

(7)

式(7)中：帶有上角標“0”的量均表示結(jié)點的值。

將應(yīng)變張量寫為Voigt向量的形式，并把式(7)中位移代入幾何方程[式(2)]，則應(yīng)變向量為

ε=Bu0

(8)

式(8)中：B為應(yīng)變矩陣。

應(yīng)力的向量形式由式(3)得到，即

σ=Dε=DBu0

(9)

式(9)中：D為彈性矩陣。

將蠕變應(yīng)變看作初應(yīng)變，忽略重力引起的體積力，根據(jù)虛功原理，可得有限元方程為

Keu0=fl+fc

(10)

Ktu0t=flt+fct

(11)

通過高斯消去法或超松弛迭代(successive over-relaxation, SOR)法即可計算t時刻的位移。下一時刻t+Δt先由式(4)和式(5)計算蠕變應(yīng)變，得到附加載荷列陣，再計算載荷向量，重新求解式(11)即得到下一時刻的位移，以此類推直至燃料棒包殼的壽期末。

1.5 包殼失效準則

為了預(yù)測包殼管服役狀態(tài)，并判斷其是否失效，需要建立包殼管的蠕變坍塌準則，擬采用以下的失效準則。

(1) 最大屈服應(yīng)力準則。當應(yīng)力達到包殼管的屈服應(yīng)力時，認為包殼失效[2]。

(2) 最大橢圓度準則和最大位移準則。當包殼管上任一點的橢圓度或位移超過限值時，認為包殼蠕變坍塌，橢圓度限值取25%[10]，最大位移的限值取包殼管最小外徑的一半[14]。

2 結(jié)果與討論

2.1 測試和驗證

為了驗證自主軟件BINE3D的正確性，擬采用相同的模型與商業(yè)軟件ABAQUS對比。分別從靜力結(jié)果和蠕變結(jié)果兩方面對自主軟件進行測試與驗證。

2.1.1 靜力測試

如圖2所示，考慮一邊長為1×1×1的立方體，在一個面受到均勻外部拉力1，材料的彈性模量為1、泊松比為0.346，理論結(jié)果、ABAQUS和本軟件結(jié)果見表2。

表2 靜力結(jié)果對比

圖2 單元示意圖

對比可看出，軟件結(jié)果與理論解和ABAQUS結(jié)果相同，BINE3D計算結(jié)果正確。

2.1.2 蠕變測試

仍采用2.1.1節(jié)的模型和參數(shù)，同時將蠕變引入，因ABAQUS軟件未提供蠕變本構(gòu)模型[式(4)]，故用經(jīng)典Norton本構(gòu)模型代替。Norton蠕變本構(gòu)的表達式為

(12)

式(12)中：常數(shù)A=0.000 713 05，n=0.8，m=-0.991 5。

BINE3D結(jié)果和商業(yè)軟件ABAQUS結(jié)果如圖3所示。

圖3 單元蠕變結(jié)果對比

可以看出，BINE3D結(jié)果和ABAQUS結(jié)果重合，具有一致性，從而驗證了BINE3D計算結(jié)果的準確性。

2.2 包殼管結(jié)構(gòu)蠕變坍塌分析

測試和驗證表明BINE3D可用于考慮蠕變效應(yīng)的三維包殼管變形分析。擬進行一真實工況下包殼管的變形分析,包殼管外徑為9.75 mm，厚度為0.57 mm，軸向高度為100 mm，彈性模量為絕對溫度的函數(shù)，即E=1.149×105-59.9T，泊松比為0.346，蠕變本構(gòu)方程為式(4)，包殼兩端位移固定，堆內(nèi)輻照參數(shù)取一典型工況下的溫度史、壓力史和中子注量史。

包殼截面為均勻圓管時，其最大位移隨著時間變化曲線見圖4。可以看出，隨著包殼管輻照時間的增加，因蠕變產(chǎn)生位移逐漸增加，并且位移的增加速度逐漸減小。到60 000 h，最大位移為0.006 81 mm，包殼結(jié)構(gòu)未失效。

圖4 包殼管最大位移隨時間變化曲線

不同初始橢圓度下三維包殼管一中部截面(高度h=42.875 mm)的橢圓度隨時間變化規(guī)律見圖5，隨著服役時間的增加，初始橢圓度較大的包殼管總橢圓度也較大。隨著輻照時間的增加，包殼橢圓度增加速度逐漸減小。輻照60 000 h時，初始橢圓度為0.1 mm的包殼最大橢圓度為3.91%，遠小于限值25%，未在服役期間出現(xiàn)坍塌。

圖5 包殼管初始橢圓度對總橢圓度的影響

當橢圓度為0.1 mm時，不同時刻包殼管一中部截面(高度h=42.875 mm)外輪廓見圖6，隨著輻照時間的增加，受外壓和持續(xù)蠕變導(dǎo)致包殼管的位移越來越大。其長軸向外變形、短軸向內(nèi)變形。

圖6 初始橢圓度為0.1 mm時，包殼管中間截面外輪廓圖

三維變形分析軟件BINE3D不需對軸向變形做出假設(shè)，也無需對包殼的變形進行修正，計算出的結(jié)果即為包殼的真實變形。隨著輻照時間的增加，長軸和短軸處軸向節(jié)點的位移見圖7。當輻照時間為60 000 h時，長軸位移呈現(xiàn)出兩端為負、中間為正的特點，這是由包殼的三維結(jié)構(gòu)導(dǎo)致的。60 000 h時，長軸最大正向位移和負向位移分別為0.002 88 mm和-0.004 mm。而基于平面應(yīng)變假設(shè)的修正方法無法精確修正出此結(jié)果。短軸位移一直為負，即短軸處向包殼內(nèi)變形，當輻照時間為60 000 h時，最大位移為-0.016 7 mm，短軸位移較長軸大。

圖7 橢圓度為0.1 mm時，包殼管長軸和短軸處軸向位置的位移隨時間變化

3 結(jié)論

基于連續(xù)介質(zhì)力學(xué)模型對燃料棒包殼進行了三維力學(xué)建模，并采用有限元方法開發(fā)出了包殼三維變形分析軟件BINE3D，可用于預(yù)測包殼的堆內(nèi)變形行為并揭示包殼的失效機理。主要工作和結(jié)論如下。

(1) 基于有限元數(shù)值方法開發(fā)了燃料棒包殼三維變形分析軟件BINE3D，可準確預(yù)測包殼三維在堆變形狀態(tài)。

(2) 包殼的變形速度(主要由蠕變導(dǎo)致)隨著輻照時間的增加逐漸減小。

(3) 具有橢圓截面的三維包殼結(jié)構(gòu)，其長軸位移有正有負，而短軸位移均為負且較長軸的變形量較大。

(4) 該軟件有望對其他復(fù)雜結(jié)構(gòu)的核燃料包殼進行分析。