壓水堆燃料棒芯塊事故后傳熱的簡化分析

齊宇博，張　偉，于　江，南金秋，趙劍剛

壓水堆燃料棒芯塊事故后傳熱的簡化分析

齊宇博，張偉，于江，南金秋，趙劍剛

（中廣核研研究院有限公司，廣東 深圳 518124）

壓水堆核電廠在嚴(yán)重事故下，堆芯換熱條件惡劣，此時包含衰變熱的堆芯非穩(wěn)態(tài)傳熱分析過程較為復(fù)雜。本文以大亞灣核電廠M310機組一回路熱段雙端斷裂為假想工況，用簡化分析方法研究事故后衰變熱的傳遞情況，獲得了含時間變量的溫度場簡化計算公式。采用ANSYS有限元分析軟件，用非簡化分析的方法計算分析，兩種方法對比驗證。研究結(jié)果顯示：在一定條件下，簡化計算的誤差相對較小，可采取簡化分析方法對事故后燃料棒芯塊的傳熱進(jìn)行計算。本文的簡化方法不僅局限于壓水堆，其他堆型的事故分析亦可借鑒。

瞬態(tài)導(dǎo)熱；集總參數(shù)法；有限元法；嚴(yán)重事故

非穩(wěn)態(tài)傳熱過程中，控制方程為幾何空間的函數(shù)，是時間尺度一階偏導(dǎo)，對較復(fù)雜工程問題，幾乎不能獲得解析解。然而，當(dāng)求解域內(nèi)熱阻較小，可近似認(rèn)為溫度場在瞬態(tài)均勻一致，此時導(dǎo)熱微分方程簡潔且易于計算，該方法稱之為集總參數(shù)法[1]。然而傳統(tǒng)的集總參數(shù)法，多數(shù)情況下回避研究對象內(nèi)部釋熱，原因是內(nèi)熱源隨時間變化，難以用確定的函數(shù)關(guān)系式表征，即使某些情況下得到表達(dá)式，控制方程依然無法解析。許多工程問題又存在內(nèi)熱源，這使得集總參數(shù)法在實際應(yīng)用過程中受到一定程度的限制。

壓水堆核電廠一回路發(fā)生大破口失水事故（LOCA）時，在緊急停堆且未啟動再循環(huán)階段，燃料棒芯塊的衰變熱近似指數(shù)衰減。本文借助origin后處理軟件，用非線性指數(shù)函數(shù)擬合該衰變關(guān)系，獲得較理想的結(jié)果。將擬合關(guān)系式引入非穩(wěn)態(tài)導(dǎo)熱微分方程中，求解并獲得分析事故后燃料芯塊瞬態(tài)導(dǎo)熱的簡化公式。同時，用ANSYS有限元軟件進(jìn)行對比驗證，在畢渥數(shù)（v）不超過0.113 75，最大誤差為13.1%。

因此，該簡化方法是可應(yīng)用于事故后燃料元件傳熱的初步分析。

1　研究的對象

壓水堆核電廠在LOCA事故后，兩臺低壓安注泵和兩臺安全殼噴淋泵投運，從安注箱和換料水箱取水。大約持續(xù)1 200 s，上述水源枯竭，安全注入系統(tǒng)（RIS）和安全殼噴淋系統(tǒng)（EAS）切換到從地坑取水，進(jìn)入長期再循環(huán)階段。

在正常運行工況下，燃料芯塊中心局部熱點溫度可達(dá)1 670 ℃[2]，為探索集總參數(shù)法在嚴(yán)重事故下應(yīng)用可行性，假設(shè)UO2熱物性參數(shù)不隨富集度、輻照、溫度等變化，忽略芯塊腫脹和變形。

1.1　物性參數(shù)

UO2芯塊物性參數(shù)[2]如表1所示。

表1　物性參數(shù)

1.2　熱邊界條件

LOCA事故后緊急停堆，假設(shè)燃料芯塊初始溫度分布均勻，取1 670 ℃；燃料棒環(huán)境溫度為90 ℃[3]。

便于后續(xù)計算和分析，燃料棒芯塊與環(huán)境間熱量交換等效成表面復(fù)合換熱系數(shù)，分別取20 W/（m2·℃）、200 W/（m2·℃）和500 W/（m2·℃）三種工況進(jìn)行參數(shù)化分析。

事故后，燃料棒芯塊衰變熱如圖1所示。

圖1　LOCA事故下芯塊體積釋熱率隨時間曲線

2　公式的推導(dǎo)

非穩(wěn)態(tài)、有內(nèi)熱源導(dǎo)熱微分方程：

忽略溫度場隨幾何空間的變化，導(dǎo)熱微分方程簡化為：

——換熱面積；

——體積；

——邊界等效換熱系數(shù)；

導(dǎo)熱微分方程可進(jìn)一步寫成：

為常數(shù)項，由初始條件進(jìn)行決定。

LOCA事故過程中，前1 200 s堆芯的衰變熱接近指數(shù)衰減，如圖2所示。

圖2　origin 擬合

非線性擬合結(jié)果：

將方程（7）帶入（6），求解：

常數(shù)項C：

根據(jù)方程（8）和（9），代入物性參數(shù)，就可得到LOCA事故后1 200 s內(nèi)，燃料芯塊瞬態(tài)導(dǎo)熱的近似解。

3　ANSYS驗證

實際中，依托計算機技術(shù)進(jìn)行數(shù)值分析，瞬態(tài)導(dǎo)熱分析常采用有限差分法（FDM）和有限元法（FEM）。首先將求解域在空間上離散，對單元體進(jìn)行高階差值或線性插值，利用計算機獲得較高計算精度和自適應(yīng)性。

瞬態(tài)熱力學(xué)分析一般方程：

式中：[]——傳導(dǎo)矩陣，包括導(dǎo)熱系數(shù)、對流系數(shù)及輻射系數(shù)和形狀系數(shù)；

[]——比熱矩陣；

使用ANSYS有限元軟件時，基于燃料芯塊對稱性，取1/4圓柱體建模，網(wǎng)格劃分如圖3所示。

圖3　網(wǎng)格劃分

3.1　表面復(fù)合換熱系數(shù)h=20 W/（m2·℃）

通常用畢渥數(shù)（v）來表征物體內(nèi)部熱阻與外部熱阻間關(guān)系，其定義：

其中：——物體長度量綱；

——物體表面?zhèn)鳠嵯禂?shù)；

芯塊的表面復(fù)合換熱系數(shù)一般不會達(dá)到該量級，但當(dāng)發(fā)生超設(shè)計基準(zhǔn)事故，例如堆芯堵塞，到達(dá)堆芯的流量不足，燃料組件上部區(qū)域處核態(tài)沸騰，換熱系數(shù)可低至該量級[3]。

簡化公式（8）和FEM計算結(jié)果如表2所示，F(xiàn)EM穩(wěn)態(tài)溫度分布云圖如圖4所示。簡化公式（8）和FEM計算結(jié)果體現(xiàn)在曲線中，如圖5所示。

表2　數(shù)據(jù)對比

注：① 指用簡化公式（8）計算結(jié)果。

② 指用FEM方法計算的結(jié)果。

③ 指用FEM方法計算芯塊表面溫度和中心溫度差值，與簡化公式（8）計算的平均溫度之比。

圖4　1 200 s FEM計算結(jié)果

結(jié)果顯示：考慮芯塊導(dǎo)熱熱阻，1 200 s后中心最高溫度為876.84 ℃，表面溫度為872.48 ℃，溫度差為3.56 ℃。相對于整體溫度，溫差非常小。

圖5　結(jié)果比較

從變化曲線看，兩種方法計算結(jié)果非常吻合，相對誤差小于0.5%。

3.2　表面復(fù)合換熱系數(shù)h=200 W/（m2·℃）

在該工況下，v=0.045 5。

簡化公式（8）和FEM計算結(jié)果如表3所示，F(xiàn)EM穩(wěn)態(tài)溫度分布云圖如圖5所示。簡化公式（8）和FEM計算結(jié)果體現(xiàn)在曲線中，如圖6所示。

表3　數(shù)據(jù)對比

注：① 指用簡化公式（8）計算結(jié)果。

② 指用FEM方法計算的結(jié)果。

③ 指用FEM方法計算芯塊表面溫度和中心溫度差值，與簡化公式（8）計算的平均溫度之比。

圖6　1 200 s FEM計算結(jié)果

結(jié)果顯示：考慮芯塊的導(dǎo)熱熱阻，在該工況下，1 200 s后中心最高溫度為166.55 ℃，表面溫度為163.22 ℃，溫度差為3.33 ℃。相對于整體溫度，該溫差并不大。此時，簡化分析方法與有限元計算結(jié)果的符合性較高。

從圖7中可見，兩種方法計算結(jié)果吻合，相對誤差小于5%，從工程角度來講是可接受的。

圖7　結(jié)果比較

3.3　表面復(fù)合換熱系數(shù)h=500 W/（m2·℃）

在該工況下，v=0.113 75。

簡化公式（8）和FEM計算結(jié)果如表4所示，F(xiàn)EM穩(wěn)態(tài)溫度分布云圖如圖8所示。簡化公式（8）和FEM計算結(jié)果體現(xiàn)在曲線中，如圖9所示。

表4　數(shù)據(jù)對比

注：① 指用簡化公式（8）計算結(jié)果。

② 指用FEM方法計算的結(jié)果。

③ 指用FEM方法計算芯塊表面溫度和中心溫度差值，與簡化公式（8）計算的平均溫度之比。

圖8　1 200 s FEM計算結(jié)果

結(jié)果顯示：考慮芯塊導(dǎo)熱熱阻，在該熱邊界條件下，1 200 s后中心最高溫度為122.61 ℃，表面溫度為119.28 ℃，溫度差為3.33 ℃。

圖9　結(jié)果比較

起始段相對誤差較大，最大相對誤差達(dá)到13.1%，60 s后兩種方法計算結(jié)果相對吻合，誤差小于5%，此時芯塊已接近穩(wěn)態(tài)。

4　結(jié)論

實際物體的幾何形狀復(fù)雜，通?？刂品匠屉y以獲得解析解，而集總參數(shù)法有效地避開幾何尺度影響，使分析過程大大簡化，當(dāng)畢渥數(shù)較小時，獲得相對精確計算結(jié)果。有國外文獻(xiàn)[3]報道，當(dāng)堆芯面臨堵塞，堆芯上部為核態(tài)沸騰區(qū)，對流換熱系數(shù)處于102量級。此時采用集總參數(shù)法進(jìn)行分析和計算是可行的，有限元計算也證明了本文所推導(dǎo)出公式有較高的計算精度。

本文將堆芯衰變功率擬合為指數(shù)形式，忽略芯塊導(dǎo)熱熱阻，求解一階非線性微分方程獲得研究芯塊導(dǎo)熱簡化公式，并采用ANSYS有限元數(shù)值解法對比驗證，當(dāng)畢渥數(shù)較小時，能夠獲得滿意的計算結(jié)果。本文的研究思路并不僅限于壓水堆事故分析，其他堆型亦可借鑒。

［1］ 楊世銘、陶文銓. 傳熱學(xué)［M］．北京：高等教育出版社，1998：66-69.

［2］ 廣東核電培訓(xùn)中心. 900 MW壓水堆核電廠系統(tǒng)與設(shè)備［M］．北京：原子能出版社，2004：50-53.

［3］ Westinghouse Electric Company. Evaluation of Long-Term Cooling Considering Particulate、Fibrous and Chemical Debris in the Recirculating Fluid，WCAP-16793-NP［R］．U.S，2011.

［4］ 黃厚誠，王秋良. 傳導(dǎo)熱問題的有限元分析［M］．北京：科學(xué)出版社，2009：1-3.

［5］ Westinghouse. Evaluation of Downstream Sump Debris Effects in Support of GSI-191，WCAP-16406-P-A［R］. U.S，2008.

Simplified Method on the Heat Conduction of PWR Fuel Rods During the Accident

QI Yubo，ZHANG Wei，YU Jiang，NAN Jinqiu，ZHAO Jiangang

（China Nuclear Power Technology Research Institute Co.，Ltd，Shenzhen of Guangdong Prov.518124，China）

During the serious accident of PWR nuclear power plants，such as a loss of coolant accident，the analysis of core decay heat transfer will be complicated。This paper is based on the hypothetical condition of a double-ended loop break of the Daya Bay M310 and uses simplified method to research the decay heat transfer，by which getting the simplified calculating formula of the temperature filed via the time variable.Meanwhile，the results of the calculation are compared with the results of ANSYS software to testify.The results have showed that the errors of calculations can be accepted when taking a simplified analytical method on a certain conditions.This simplified analytical method not only can be used for PWR，but is also available for the other types of reactors.

Transient conduction；Lumped parameter method；Finite element method；Severe accident

TLLT331

AA

0258-0918（2022）01-0053-06

2020-09-24

國家科技支撐計劃（2011BAA06B01）資助

齊宇博（1986—），男，陜西定邊人，工程師，學(xué)士，現(xiàn)主要從事核電、氫能方面研究