考慮各向異性蠕變的鋯包殼鼓脹行為數(shù)值模擬方法研究

李 偉,李小雨,段倩妮,王皓坤,武俊梅,劉仕超

(1.西安交通大學(xué) 復(fù)雜服役環(huán)境重大裝備結(jié)構(gòu)強(qiáng)度與壽命全國(guó)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710049;2.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室,四川 成都 610213)

在反應(yīng)堆運(yùn)行瞬態(tài)或事故工況的高溫環(huán)境中,Zr合金包殼的強(qiáng)度急劇降低,蠕變速率加快,在裂變氣體造成的管內(nèi)、外壓差作用下將出現(xiàn)鼓脹和爆破失效。包殼鼓脹和爆破對(duì)反應(yīng)堆安全帶來(lái)一系列不利影響,特別是對(duì)于高燃耗燃料,包殼鼓脹使燃料芯塊失去徑向約束,燃料碎塊容易重定位至鼓脹區(qū)域而使包殼局部熱集中,鼓脹區(qū)域因此被加速氧化。國(guó)內(nèi)和國(guó)際上針對(duì)輕水堆燃料包殼的失水事故(LOCA)行為開(kāi)展了大量實(shí)驗(yàn)研究[1-5]?？紤]到包殼鼓脹和爆破失效行為的復(fù)雜性,通過(guò)發(fā)展先進(jìn)、可靠的計(jì)算工具進(jìn)行預(yù)測(cè)分析有助于加強(qiáng)對(duì)該科學(xué)問(wèn)題的認(rèn)識(shí),理解包殼性能對(duì)堆芯安全的影響,有利于開(kāi)展事故的最佳估算分析。美國(guó)從2008年開(kāi)始開(kāi)發(fā)了有限元燃料性能分析程序BISON,可用于LOCA下的燃料元件行為分析[6],具備二維軸對(duì)稱(chēng)和三維建模能力。但是,BISON未考慮α相各向異性蠕變。韓國(guó)原子能研究院(KAERI)研究人員采用簡(jiǎn)化的軸對(duì)稱(chēng)解析解對(duì)包殼的LOCA高溫?zé)?力學(xué)行為進(jìn)行了分析[7-8],表明由于α相的各向異性蠕變行為,軸向單軸拉伸條件下獲得的蠕變模型與有限元方法中采用的等效應(yīng)力蠕變模型存在較大差別。

針對(duì)α相區(qū)間的Zr包殼各向異性蠕變問(wèn)題,本文基于Hill準(zhǔn)則推導(dǎo)得到應(yīng)力更新算法和一致切線剛度更新算法?；谟邢拊浖嗀BAQUS的材料子程序UMAT,通過(guò)加入相轉(zhuǎn)變、各向異性的高溫蠕變和失效限值經(jīng)驗(yàn)?zāi)Ｐ?初步獲得精細(xì)化的鼓脹行為分析工具,通過(guò)解析解和PUZRY系列爆破實(shí)驗(yàn)對(duì)該分析工具的合理性和可靠性進(jìn)行了確認(rèn)和驗(yàn)證,可為后續(xù)高燃耗燃料元件及耐事故燃料包殼(如涂層Zr和FeCrAl包殼)的事故失效行為研究提供參考。

1 材料模型

本研究采用的Zr-4合金導(dǎo)熱系數(shù)、比熱容、楊氏模量、泊松比、熱膨脹系數(shù)等基本熱-力學(xué)物性模型均取自MATPRO手冊(cè)[9]。相變模型來(lái)自文獻(xiàn)[10],該模型描述了在恒定加熱或冷卻速率下β相的體積份額。此外,本文采用兩種經(jīng)驗(yàn)?zāi)Ｐ蛠?lái)判斷是否達(dá)到爆破,即:FRAPTRAN程序應(yīng)變限值[11]和Rosinger應(yīng)力限值[12]。在計(jì)算過(guò)程中,當(dāng)包殼管中的等效蠕變應(yīng)變或環(huán)向真應(yīng)力達(dá)到任何一個(gè)限值時(shí)認(rèn)為包殼管發(fā)生爆破并終止計(jì)算。Rosinger爆破應(yīng)力模型包括3個(gè),分別為上限、平均和下限爆破應(yīng)力模型。

表1 公開(kāi)文獻(xiàn)中的高溫蠕變模型

表1中的蠕變速率參數(shù)A、n、Q與β相的體積份額有關(guān)。已知純?chǔ)料唷?0% α相和50% β混合相以及純?chǔ)孪嗟娜渥兯俾蕝?shù),本文采用傳統(tǒng)的加權(quán)混合計(jì)算方法[4],依據(jù)β相的體積份額對(duì)表1中的參數(shù)Q和n進(jìn)行線性插值,對(duì)系數(shù)A的自然對(duì)數(shù)lnA進(jìn)行線性差值,計(jì)算得到對(duì)應(yīng)的蠕變速率參數(shù)。

2 蠕變隱式積分算法

盡管Zr包殼的彈性變形也是各向異性[15],但各向異性程度較小。為簡(jiǎn)化計(jì)算,本文假設(shè)Zr包殼的彈性力學(xué)性質(zhì)為各向同性,因此只需兩個(gè)獨(dú)立的彈性常數(shù),如通常采用楊氏模量和泊松比。對(duì)于Zr-4合金,本文采用MATPRO手冊(cè)[9]中的楊氏模量和泊松比模型:

(1)

K1=(6.61×105+5.912×102T)·

(ηZr-1.2×10-3)

(2)

K2=-2.6×104CWZr

(3)

(4)

νZr=0.330 3+8.375×10-5(T-273.15)

(5)

式中:EZr為楊氏模型,MPa;T為溫度,K;CWZr為冷加工量;ηZr為含氧量,用質(zhì)量分?jǐn)?shù)表示;Φ為中子注量,m-2;νZr為泊松比。

本文采用Hill準(zhǔn)則[16]描述α相Zr合金的高溫蠕變變形為各向異性。定義一個(gè)屈服函數(shù)φ為:

(6)

式中:σ為柯西應(yīng)力張量,為了方便數(shù)值算法的編程實(shí)現(xiàn),本文采用列陣形式,對(duì)于三維單元,σ=[σ11,σ22,σ33,σ12,σ13,σ23]T;P為對(duì)稱(chēng)矩陣,對(duì)于三維幾何,其形式為:

(7)

式中,F、G、H、L、M、N為材料的各向異性常數(shù),針對(duì)特定材料,需通過(guò)實(shí)驗(yàn)加以確定。

Hill等效應(yīng)力與屈服函數(shù)的關(guān)聯(lián)如下:

(8)

可以看出,當(dāng)各向異性常數(shù)設(shè)置為F=G=H=0.5時(shí),式(8)即為Von Mises應(yīng)力。根據(jù)相關(guān)聯(lián)流動(dòng)法則,蠕變應(yīng)變分量與等效蠕變應(yīng)變的關(guān)系為:

(9)

基于有限元方法的蠕變隱式積分算法包含兩個(gè)步驟:局部迭代和一致切線剛度。局部迭代用于實(shí)現(xiàn)應(yīng)力的更新,需滿(mǎn)足的等式如下:

(10)

(11)

為了進(jìn)一步展開(kāi)式(11),將時(shí)間步末的應(yīng)力張量寫(xiě)為:

(12)

其中:

(13)

(14)

(15)

σtr=σt+DeΔε

(16)

式中:I為單位矩陣;De為彈性矩陣;Δε為總應(yīng)變張量的增量;Δεcr為蠕變應(yīng)變張量的增量;σtr為試應(yīng)力張量;σt為上個(gè)時(shí)間步末的應(yīng)力張量,其余變量為運(yùn)算所需的中間變量。

結(jié)合式(6)可求得應(yīng)變張量對(duì)等效蠕變應(yīng)變?cè)隽康钠珜?dǎo)數(shù)為:

(17)

當(dāng)各向異性常數(shù)為F=G=H=0.5、L=M=N=1.5時(shí),可以證明式(17)等于-3μ,其中μ為剪切模量,這與各向同性Von Mises的情況相符。當(dāng)上述局部迭代達(dá)到收斂后,根據(jù)式(13),應(yīng)力也得以更新。

對(duì)于一致切線剛度矩陣,其定義為:

(18)

式中:J為一致切線剛度矩陣;Δσ為應(yīng)力張量的增量。

對(duì)式(13)兩邊求偏導(dǎo),并考慮式(9)得到:

(19)

(20)

經(jīng)過(guò)進(jìn)一步整理,可得到一致切線剛度矩陣為:

(21)

其中:

(22)

(23)

(24)

3 方法驗(yàn)證

3.1 厚壁圓管蠕變解析解

a——圓管內(nèi)表面等效蠕變應(yīng)變隨時(shí)間的變化;b——圓管徑向上的應(yīng)力分布,圓管內(nèi)表面等效蠕變應(yīng)變?yōu)?.66

3.2 PUZRY實(shí)驗(yàn)

采用PUZRY系列鼓脹爆破實(shí)驗(yàn)[18]進(jìn)行驗(yàn)證。包殼管長(zhǎng)50 mm,兩端各有長(zhǎng)5 mm的端塞區(qū),內(nèi)徑9.3 mm,外徑10.75 mm。材料為未輻照、未氧化的Zr-4。PUZRY實(shí)驗(yàn)在氦氣氛圍中進(jìn)行,包含31個(gè)工況,溫度范圍為698～1 201 ℃,涵蓋了α相、α+β混合相和β相。本文模擬了其中的28個(gè)工況,未模擬爆破時(shí)間超過(guò)4 000 s的工況。為減少計(jì)算時(shí)間,利用了軸對(duì)稱(chēng)幾何。對(duì)于α相,各向異性常數(shù)取F=0.95、G=0.304、H=0.24、L=M=N=1.5[19]。僅對(duì)α相區(qū)間的蠕變施加各向異性蠕變算法,對(duì)于α+β混合相和β相,認(rèn)為其為各向同性蠕變。

圖2示出PUZRY實(shí)驗(yàn)數(shù)據(jù)與FRAPTRAN爆破應(yīng)變模型和Rosinger爆破應(yīng)力模型的對(duì)比,其中爆破應(yīng)力實(shí)驗(yàn)數(shù)據(jù)根據(jù)內(nèi)壓和環(huán)向應(yīng)變實(shí)驗(yàn)數(shù)據(jù)計(jì)算所得[20]。由圖2可見(jiàn):FRAPTRAN爆破應(yīng)變模型在α+β混合相和β相區(qū)間與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比結(jié)果不佳,而在α相區(qū)間兩者的符合程度稍好;Rosinger平均爆破應(yīng)力模型與α+β混合相區(qū)間的實(shí)驗(yàn)數(shù)據(jù)符合較好;Rosinger下限爆破應(yīng)力模型與β相區(qū)間的實(shí)驗(yàn)數(shù)據(jù)符合較好。因此,本文在這3個(gè)區(qū)域分別采用了FRAPTRAN爆破應(yīng)變模型、Rosinger平均爆破應(yīng)力模型和Rosinger下限爆破應(yīng)力模型來(lái)判斷實(shí)驗(yàn)樣品是否發(fā)生了爆破。

a——爆破應(yīng)變模型;b——爆破應(yīng)力模型

計(jì)算得到的爆破時(shí)間與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖3所示。除了工況6和19,程序預(yù)測(cè)的爆破時(shí)間與實(shí)驗(yàn)數(shù)據(jù)均符合較好。這兩個(gè)工況的溫度分別為約950 ℃和900 ℃,位于α+β混合相區(qū)間,表明本文根據(jù)文獻(xiàn)中的簡(jiǎn)單加權(quán)混合方法計(jì)算該區(qū)間的蠕變應(yīng)變率存在較大偏差,而應(yīng)開(kāi)發(fā)更為機(jī)理性的模型。由圖3可見(jiàn),相對(duì)于各向同性蠕變算法,各向異性蠕變算法能明顯提高對(duì)α相區(qū)間蠕變速率預(yù)測(cè)的準(zhǔn)確度,這是因?yàn)槿渥兯俾誓Ｐ屯ǔ；趩屋S應(yīng)力實(shí)驗(yàn)數(shù)據(jù)獲得,在將單軸應(yīng)力用等效應(yīng)力代替時(shí),各向異性情況下的蠕變速率模型需要根據(jù)各向異性常數(shù)進(jìn)行調(diào)整,從而改變了蠕變速率。此外,各向異性情況下,等效應(yīng)力也與各向同性情況下采用的Von Mises應(yīng)力值不同。圖4示出不同溫度下計(jì)算得到爆破時(shí)間與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比,直觀顯示了該算法對(duì)α相區(qū)間各向異性蠕變速率的適用性。圖4結(jié)果還表明:對(duì)于PUZRY實(shí)驗(yàn)條件,Rosinger蠕變模型對(duì)于β相區(qū)間的適用性較好。

a——各向異性蠕變算法;b——各向同性蠕變算法

a——各向異性蠕變算法;b——各向同性蠕變算法

爆破時(shí)外表面環(huán)向應(yīng)變與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖5所示。本文未采用單一的爆破準(zhǔn)則,而是針對(duì)不同的相區(qū)間采用不同的爆破準(zhǔn)則,使得計(jì)算結(jié)果在總體上與實(shí)驗(yàn)數(shù)據(jù)符合較好。但在β相區(qū)間,對(duì)于加壓速率大的工況(7.1×10-3MPa/s),程序預(yù)測(cè)的爆破環(huán)向應(yīng)變較為明顯低于實(shí)驗(yàn)數(shù)據(jù),對(duì)于加壓速率小的工況(6.4×10-4MPa/s),計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)符合較好。這可能是由于當(dāng)加壓速率較大時(shí),塑性變形的作用不可忽視,但現(xiàn)階段僅考慮了蠕變,導(dǎo)致應(yīng)力的預(yù)測(cè)值過(guò)大而過(guò)早爆破。

圖5 爆破應(yīng)變計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)對(duì)比

圖6示出包殼軸向變形(實(shí)驗(yàn)后的包殼長(zhǎng)度/包殼初始長(zhǎng)度)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比。由于包殼鼓脹造成的環(huán)向伸長(zhǎng),在總體積保持不變的情況下,軸向?qū)l(fā)生一定程度的收縮。在α+β混合相區(qū)間和β相區(qū)間,程序預(yù)測(cè)的軸向變形與實(shí)驗(yàn)數(shù)據(jù)符合較好(注意在這兩個(gè)相區(qū)間均為各向同性蠕變),而在α相區(qū)間,包殼鼓脹最為明顯,但采用各向同性蠕變算法預(yù)測(cè)的軸向收縮量明顯偏低。在各向同性和各向異性蠕變?cè)讦料鄥^(qū)間的爆破環(huán)向應(yīng)變相同(均采用了FRAPTRAN爆破應(yīng)變準(zhǔn)則)、而軸向變形差異明顯的情況下,表明采用各向同性蠕變算法計(jì)算α相的蠕變時(shí),包殼厚度必然比采用各向異性蠕變算法時(shí)小,這與數(shù)值模擬結(jié)果相符。

a——各向同性蠕變算法;b——各向異性蠕變算法

圖7示出包殼軸向中間位置處的環(huán)向應(yīng)變隨時(shí)間的典型變化。圖7中,PUZRY后的編號(hào)表示PUZRY實(shí)驗(yàn)的工況號(hào),例如PUZRY-26表示第26實(shí)驗(yàn)工況。由圖7可見(jiàn),對(duì)于α相區(qū)間,各向異性蠕變算法比各向同性蠕變算法預(yù)測(cè)的包殼環(huán)向應(yīng)變變化更慢,在這種情況下,預(yù)測(cè)的準(zhǔn)確度更依賴(lài)于可靠的爆破模型。但對(duì)于α+β混合相區(qū)間,蠕變模型本身帶來(lái)的蠕變速率過(guò)大,造成計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)差異明顯,因此有必要開(kāi)發(fā)更機(jī)理性的蠕變模型。對(duì)于β相區(qū)間,預(yù)測(cè)結(jié)果的不確定性似乎還與壓力加速率大小有關(guān)。圖8示出工況PUZRY-30(α相區(qū)間)的計(jì)算結(jié)果及與實(shí)驗(yàn)[21]的對(duì)比。由圖8可見(jiàn),各向異性蠕變算法預(yù)測(cè)的包殼發(fā)生鼓脹的軸向區(qū)間比各向同性蠕變算法更大(即后者預(yù)測(cè)的鼓脹更為集中在包殼的中間位置),各向異性蠕變算法預(yù)測(cè)的鼓脹平均環(huán)向應(yīng)變比各向同性蠕變算法預(yù)測(cè)的值更大。

a——PUZRY-26;b——PUZRY-30;c——PUZRY-18;d——PUZRY-12,PUZRY-1

a——各向同性蠕變算法;b——各向異性蠕變算法

4 結(jié)論

在高溫環(huán)境和內(nèi)外壓差條件下,Zr包殼不可避免的存在鼓脹爆破失效。對(duì)此現(xiàn)象進(jìn)行高保真數(shù)值模擬研究在高燃耗背景下顯得越來(lái)越重要,因?yàn)榘鼩さ墓拿洷剖c高燃耗下的反應(yīng)堆安全密切相關(guān)。本文發(fā)展了α相Zr包殼的各向異性蠕變的隱式數(shù)值積分算法,且在有限元方法中加以實(shí)現(xiàn),通過(guò)厚壁圓管蠕變問(wèn)題的理論解析解進(jìn)行了正確性確認(rèn),結(jié)合FRAPTRAN爆破應(yīng)變限值模型和Rosinger爆破應(yīng)力限值模型,模擬了PUZRY鼓脹爆破實(shí)驗(yàn),開(kāi)展了模型和算法的合理性驗(yàn)證。結(jié)果表明:使用各向異性蠕變算法比各向同性蠕變算法能更好地預(yù)測(cè)α相區(qū)間的鼓脹;對(duì)于α+β混合相區(qū)間,目前常用的加權(quán)混合方法計(jì)算蠕變速率還存在較大誤差;對(duì)于β相區(qū)間,目前僅考慮了蠕變,因此對(duì)加壓速率較低的情況適用性較好,但當(dāng)加壓速率更高時(shí),未包含塑性應(yīng)變可能導(dǎo)致應(yīng)力預(yù)測(cè)值偏大而過(guò)早地判斷發(fā)生爆破。此外,與溫度等相關(guān)的各向異性常數(shù)、包含塑性損傷等將在后續(xù)的工作中加以考慮。