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

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

(1.西安交通大學 復雜服役環境重大裝備結構強度與壽命全國重點實驗室,陜西 西安 710049;2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室,四川 成都 610213)

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

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

1 材料模型

本研究采用的Zr-4合金導熱系數、比熱容、楊氏模量、泊松比、熱膨脹系數等基本熱-力學物性模型均取自MATPRO手冊[9]。相變模型來自文獻[10],該模型描述了在恒定加熱或冷卻速率下β相的體積份額。此外,本文采用兩種經驗模型來判斷是否達到爆破,即:FRAPTRAN程序應變限值[11]和Rosinger應力限值[12]。在計算過程中,當包殼管中的等效蠕變應變或環向真應力達到任何一個限值時認為包殼管發生爆破并終止計算。Rosinger爆破應力模型包括3個,分別為上限、平均和下限爆破應力模型。

表1 公開文獻中的高溫蠕變模型

表1中的蠕變速率參數A、n、Q與β相的體積份額有關。已知純α相、50% α相和50% β混合相以及純β相的蠕變速率參數,本文采用傳統的加權混合計算方法[4],依據β相的體積份額對表1中的參數Q和n進行線性插值,對系數A的自然對數lnA進行線性差值,計算得到對應的蠕變速率參數。

2 蠕變隱式積分算法

盡管Zr包殼的彈性變形也是各向異性[15],但各向異性程度較小。為簡化計算,本文假設Zr包殼的彈性力學性質為各向同性,因此只需兩個獨立的彈性常數,如通常采用楊氏模量和泊松比。對于Zr-4合金,本文采用MATPRO手冊[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為含氧量,用質量分數表示;Φ為中子注量,m-2;νZr為泊松比。

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

(6)

式中:σ為柯西應力張量,為了方便數值算法的編程實現,本文采用列陣形式,對于三維單元,σ=[σ11,σ22,σ33,σ12,σ13,σ23]T;P為對稱矩陣,對于三維幾何,其形式為:

(7)

式中,F、G、H、L、M、N為材料的各向異性常數,針對特定材料,需通過實驗加以確定。

Hill等效應力與屈服函數的關聯如下:

(8)

可以看出,當各向異性常數設置為F=G=H=0.5時,式(8)即為Von Mises應力。根據相關聯流動法則,蠕變應變分量與等效蠕變應變的關系為:

(9)

基于有限元方法的蠕變隱式積分算法包含兩個步驟:局部迭代和一致切線剛度。局部迭代用于實現應力的更新,需滿足的等式如下:

(10)

(11)

為了進一步展開式(11),將時間步末的應力張量寫為:

(12)

其中:

(13)

(14)

(15)

σtr=σt+DeΔε

(16)

式中:I為單位矩陣;De為彈性矩陣;Δε為總應變張量的增量;Δεcr為蠕變應變張量的增量;σtr為試應力張量;σt為上個時間步末的應力張量,其余變量為運算所需的中間變量。

結合式(6)可求得應變張量對等效蠕變應變增量的偏導數為:

(17)

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

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

(18)

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

對式(13)兩邊求偏導,并考慮式(9)得到:

(19)

(20)

經過進一步整理,可得到一致切線剛度矩陣為:

(21)

其中:

(22)

(23)

(24)

3 方法驗證

3.1 厚壁圓管蠕變解析解

a——圓管內表面等效蠕變應變隨時間的變化;b——圓管徑向上的應力分布,圓管內表面等效蠕變應變為0.66

3.2 PUZRY實驗

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

圖2示出PUZRY實驗數據與FRAPTRAN爆破應變模型和Rosinger爆破應力模型的對比,其中爆破應力實驗數據根據內壓和環向應變實驗數據計算所得[20]。由圖2可見:FRAPTRAN爆破應變模型在α+β混合相和β相區間與實驗數據的對比結果不佳,而在α相區間兩者的符合程度稍好;Rosinger平均爆破應力模型與α+β混合相區間的實驗數據符合較好;Rosinger下限爆破應力模型與β相區間的實驗數據符合較好。因此,本文在這3個區域分別采用了FRAPTRAN爆破應變模型、Rosinger平均爆破應力模型和Rosinger下限爆破應力模型來判斷實驗樣品是否發生了爆破。

a——爆破應變模型;b——爆破應力模型

計算得到的爆破時間與實驗數據的對比如圖3所示。除了工況6和19,程序預測的爆破時間與實驗數據均符合較好。這兩個工況的溫度分別為約950 ℃和900 ℃,位于α+β混合相區間,表明本文根據文獻中的簡單加權混合方法計算該區間的蠕變應變率存在較大偏差,而應開發更為機理性的模型。由圖3可見,相對于各向同性蠕變算法,各向異性蠕變算法能明顯提高對α相區間蠕變速率預測的準確度,這是因為蠕變速率模型通?；趩屋S應力實驗數據獲得,在將單軸應力用等效應力代替時,各向異性情況下的蠕變速率模型需要根據各向異性常數進行調整,從而改變了蠕變速率。此外,各向異性情況下,等效應力也與各向同性情況下采用的Von Mises應力值不同。圖4示出不同溫度下計算得到爆破時間與實驗數據的對比,直觀顯示了該算法對α相區間各向異性蠕變速率的適用性。圖4結果還表明:對于PUZRY實驗條件,Rosinger蠕變模型對于β相區間的適用性較好。

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

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

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

圖5 爆破應變計算值與實驗數據對比

圖6示出包殼軸向變形(實驗后的包殼長度/包殼初始長度)與實驗數據對比。由于包殼鼓脹造成的環向伸長,在總體積保持不變的情況下,軸向將發生一定程度的收縮。在α+β混合相區間和β相區間,程序預測的軸向變形與實驗數據符合較好(注意在這兩個相區間均為各向同性蠕變),而在α相區間,包殼鼓脹最為明顯,但采用各向同性蠕變算法預測的軸向收縮量明顯偏低。在各向同性和各向異性蠕變在α相區間的爆破環向應變相同(均采用了FRAPTRAN爆破應變準則)、而軸向變形差異明顯的情況下,表明采用各向同性蠕變算法計算α相的蠕變時,包殼厚度必然比采用各向異性蠕變算法時小,這與數值模擬結果相符。

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

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

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

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

4 結論

在高溫環境和內外壓差條件下,Zr包殼不可避免的存在鼓脹爆破失效。對此現象進行高保真數值模擬研究在高燃耗背景下顯得越來越重要,因為包殼的鼓脹爆破失效與高燃耗下的反應堆安全密切相關。本文發展了α相Zr包殼的各向異性蠕變的隱式數值積分算法,且在有限元方法中加以實現,通過厚壁圓管蠕變問題的理論解析解進行了正確性確認,結合FRAPTRAN爆破應變限值模型和Rosinger爆破應力限值模型,模擬了PUZRY鼓脹爆破實驗,開展了模型和算法的合理性驗證。結果表明:使用各向異性蠕變算法比各向同性蠕變算法能更好地預測α相區間的鼓脹;對于α+β混合相區間,目前常用的加權混合方法計算蠕變速率還存在較大誤差;對于β相區間,目前僅考慮了蠕變,因此對加壓速率較低的情況適用性較好,但當加壓速率更高時,未包含塑性應變可能導致應力預測值偏大而過早地判斷發生爆破。此外,與溫度等相關的各向異性常數、包含塑性損傷等將在后續的工作中加以考慮。