張小英 姚婷婷 李志威 黃 凱

1（華南理工大學 電力學院 廣州 510640）

2（中科華核電技術(shù)研究院 深圳 518026）

堆芯熔融物對壓力容器壁面燒蝕過程的數(shù)值模擬

張小英1姚婷婷1李志威1黃 凱2

1（華南理工大學 電力學院 廣州 510640）

2（中科華核電技術(shù)研究院 深圳 518026）

研究堆芯熔融物對壓力容器壁面的動態(tài)燒蝕，對于反應堆冷卻劑嚴重喪失事故(Loss of coolant accident, LOCA)后果的預測以及緩解方案的設(shè)計具有重要意義。本文以AP600為研究對象，在假設(shè)冷卻劑全部喪失事故工況下，采用堆芯熔融物兩層結(jié)構(gòu)模型，計算熔池對壁面的加熱；建立壓力容器壁面的非穩(wěn)態(tài)二維傳熱模型，并考慮安全殼水池對壓力容器外側(cè)的冷卻，采用移動邊界模型模擬燒蝕引起壁面局部厚度變??；計算了堆芯熔融物坍塌后15000s范圍內(nèi)，壓力容器下封頭壁面溫度和厚度的變化。

反應堆，冷卻劑喪失，堆芯熔融物，壓力容器，壁面燒蝕

在反應堆發(fā)生冷卻劑嚴重喪失事故(Loss of coolant accident, LOCA)[1]后，堆芯由于長時間得不到冷卻而發(fā)生熔化，堆芯熔融物跌落到壓力容器底部，使底部結(jié)構(gòu)受到持續(xù)強烈的加熱而燒蝕，破壞壓力容器底部結(jié)構(gòu)的完整性。因此，研究堆芯熔融物對壓力容器下封頭材料的燒蝕問題，對于預測反應堆冷卻劑嚴重喪失事故的發(fā)生以及研究相應的緩解方案都有重要的參考意義。對此， 采用熔池分層結(jié)構(gòu)研究了下封頭壁面的熱負荷[2]，并假定材料厚度不變情況下，采用二維導熱模型計算壁面溫度；基于熔池兩層結(jié)構(gòu)模型，計算了對下封頭壁面的對流與輻射熱流[3]，用一體化嚴重事故計算程序計算了小破口和大破口事故下，下封頭壁面的熱流密度和熔池結(jié)構(gòu)[4]。

熾熱的堆芯熔融物沉積在壓力容器下封頭，強烈加熱壁面材料，引起壁面局部燒蝕，危及壓力容器的完整性。對于堆芯熔融物加熱作用下，壓力容器壁面的溫度上升及燒蝕過程還有待研究。為此，本文建立了熔池下封頭壁面?zhèn)鳠崮Ｐ?，其中熔池采用兩層結(jié)構(gòu)模型，壁面?zhèn)鳠岵捎梅欠€(wěn)態(tài)二維傳熱模型，熔池對壁面加熱考慮對流和輻射換熱，非穩(wěn)態(tài)計算中壁面燒蝕引起局部厚度的變化采用移動邊界模型。在假設(shè)壓力容器沉浸于安全殼水池中，外壁溫度接近安全殼環(huán)境壓力下的飽和水溫度情況下，計算了AP600反應堆[5–7]堆芯熔融物坍塌后15000s范圍內(nèi)，壁面溫度及厚度的變化。

1 熔池兩層結(jié)構(gòu)模型

反應堆堆芯熔融后，形成鈾氧化物和金屬熔融物，由于鈾氧化物密度大，其熔融物沉積在熔池的下部，金屬熔融物的密度較輕，浮升在熔池的上部，故采用圖1所示的兩層結(jié)構(gòu)模型[8]，模擬堆芯熔融物在壓力容器下封頭的滯留狀態(tài)。鈾氧化物熔池周邊由于與溫度較低的壓力容器壁面和金屬熔融層接觸，形成薄薄一層凝固殼層(圖2)。根據(jù)堆芯熔融物的體積，由于下封頭形狀近似為球形，從球冠的體積可依據(jù)式(1)計算熔融層高度，得到鈾氧化物熔池的高度HU和金屬熔融層的高度HM，及鈾氧化物熔池和金屬熔融層上表面的極角Uθ和Mθ。

圖1 堆芯熔融物在壓力容器下封頭形成的兩層結(jié)構(gòu)熔池Fig.1 Two-layer structure of melt core in the lower head of pressure vessel.

圖2 壁面網(wǎng)格離散化示意圖Fig.2 Grid discretization of the pressure vessel wall.

2 壓力容器壁面?zhèn)鳠岬挠嬎隳Ｐ?/h2>

2.1 微分方程及離散

為計算堆芯熔融物對壓力容器壁面的燒蝕過程，考慮到壓力容器下封頭的球面形狀，幾何上具有軸對稱性，因此對壓力容器下封頭壁面的傳熱可采用沿壁面厚度和下封頭圓周角的二維球面導熱模型，見圖3所示，AB段為熔池以上的部分，BC段為金屬熔池部分，CD段為鈾氧化物熔池部分。壁面導熱的微分方程為：

對圖3所示的下封頭壁面沿半徑和圓周角方向等間隔劃分控制體，采用控制容積積分法得到離散方程為[9]：

其中：

在熔池對壓力容器壁面的燒蝕過程中，壓力容器壁面局部因燒蝕而變薄，故壁面?zhèn)鳠岬挠嬎銋^(qū)域?qū)l(fā)生變化[10]。為模擬壁面燒蝕過程中計算區(qū)域的變化，本文采用移動邊界模型。對每一個新的計算時刻，首先判斷前一時刻的壁面溫度，將超過材料熔點的那些節(jié)點從計算域中去掉，采用未熔的節(jié)點構(gòu)造新的計算域，如圖4所示。

圖3 壓力容器壁面?zhèn)鳠嵊嬎銕缀文Ｐ?(a) 徑向，(b) 局部Fig.3 Geometry model of heat analysis of pressure vessel wall. (a) Radial, (b) Part

圖4 移動邊界模型示意圖Fig.4 Moving boundary model scheme.

2.2 計算的邊界條件

對于邊界控制體的離散方程，采用補充虛擬邊界控制體方法，并令虛擬邊界控制體的溫度始終等于相鄰邊界控制體的溫度，由此邊界控制體與虛擬控制體之間雖有導熱項卻沒有導熱熱流，式(3)對邊界控制體可適用。在壁面內(nèi)側(cè)，需計算壁面與熔池的對流與輻射換熱。壁面外側(cè)浸泡在安全殼水池中，近壁區(qū)的水受到壁面加熱達到沸騰，因此可令外壁溫度接近安全殼壓力下水的飽和溫度。

圖3所示的壓力容器內(nèi)側(cè)邊界各部分的受熱不同，需分別計算。AB段的內(nèi)邊界熱流需計算金屬熔池上表面對壁面的輻射換熱；BC段的內(nèi)邊界熱流需計算金屬熔池對壁面的對流換熱凝固殼層的釋熱；CD段的內(nèi)邊界熱流需計算鈾氧化物熔池對壁面的對流換熱。

為計算AB段內(nèi)邊界的輻射熱流qW,r，只考慮金屬熔融層上表面對壁面的輻射，忽略壓力容器的自輻射反射，根據(jù)Lambert定律計算金屬熔融層對壓力容器內(nèi)表面的凈輻射熱流：

BC段內(nèi)邊界的熱流為金屬熔融層對壁面的對流熱流MWq，計算為：

式中，BT是金屬熔融層主體溫度，K，需根據(jù)金屬熔融層的換熱計算。

CD段內(nèi)邊界的熱流包括鈾氧化物熔池對壁面的對流換熱qU,dn(θ)及凝固殼層的釋熱：

根據(jù)文獻[2]，鈾氧化物熔池由于溫度分層形成上對流層和下對流層，上對流層與金屬熔融層對流換熱，下對流層與壁面對流換熱，換熱熱流計算為：式中，U,dnq為下對流層與所接觸壁面的平均換熱熱流。

式中，'Ra是內(nèi)部Rayleigh數(shù)。

2.3 金屬熔融層主體溫度TB的計算

在前面計算金屬熔池對壓力容器內(nèi)壁的對流與輻射換熱當中，均用到金屬熔池的主體溫度TB。為計算TB，對金屬熔融層列熱平衡方程，單位時間鈾氧化物熔池通過對流傳給上凝固殼層，并進一步傳給金屬熔池的熱量MUq，等于金屬熔池向壓力容器壁面的輻射換熱量,Wrq和對流換熱量MWq之和，即：

式中，AM是熔池的橫截面積；AW是與金屬熔池接觸的壓力容器內(nèi)壁面積；FM–W是金屬熔池對壓力容器AB段內(nèi)壁的輻射角系數(shù)。鈾氧化物熔池對金屬熔融層的對流換熱計算采用Globe Dropkin的關(guān)系式計算[11]：

將式(5)、(6)、(12)代入式(11)，TW,in采用前一計算時刻的數(shù)值，求解可得到金屬熔池的主體溫度TB。

2.4 凝固殼層厚度的計算

凝固殼層厚度是熔池結(jié)構(gòu)的重要參數(shù)，為計算凝固殼層的厚度，采用一維平板導熱模型，對凝固殼層和壓力容器壁面沿厚度方向列熱平衡方程，即凝固殼層的徑向?qū)釤崃?壓力容器的徑向?qū)釤崃?熔池對壁面的對流熱流+堆芯的衰變熱流，由此求解得到凝固殼層的厚度。

3 計算結(jié)果與分析

在計算壓力容器下封頭壁面溫度時，取計算段高度為下封頭至金屬熔融層以上0.5m高度的壁面，考慮幾何的對稱性，取一個圓周方向的壁面加以計算，有關(guān)參數(shù)如下：鈾氧化物密度ρU=8450kg·m?3；鈾氧化物熔點TU=2973K；鈾氧化物熱導率λU=5.3W·(m·K)?1；鈾氧化物體積膨脹系數(shù)βU=1.05×10?4K?1；鈾氧化物動力粘性系數(shù)鈾氧化物發(fā)射率εU=0.8；AP600下封頭半徑R=2m；鈾氧化物熔池體積VU=8.6 m3；金屬熔池體積VM=3.14m3；金屬熔融物密度ρM=6890kg·m?3；金屬熔融物熔點TM=1600K；金屬熔融物熱導率λM=25.5W·(m·K)?1；金屬熔融物體積膨脹系數(shù)βM=1.1×10?4K?1；金屬熔融物動力粘性系數(shù)金屬熔融物發(fā)射率εM=0.45；AP600壓力容器壁厚δV=0.015m；AP600停堆4 h后，單位體積釋熱量壓力容器壁面初始溫度322K；安全殼工作壓力0.1MPa。

本文計算堆芯熔融物對下封頭壁面內(nèi)側(cè)的加熱熱流如圖5(a)所示，相應地文獻[2]計算的結(jié)果如圖5(b)所示，文獻給出了將金屬熔融層以上壁面視為絕熱（實線）和考慮金屬熔融層以上壁面?zhèn)鳠岬膬山M結(jié)果（虛線）?？梢钥吹剑?1) 在下封頭與鈾氧化物熔池接觸的部分，熔池對壁面的加熱熱流隨傾角而上升；(2) 在下封頭與金屬熔融層接觸的部分，加熱熱流下降，因為在金屬熔融物以上部分的壁面主要的加熱來自金屬熔融層表面的輻射，熱流比金屬熔融層以下壁面受到的對流加熱小得多，導致金屬熔融層接觸部分的壁面加熱熱流也有所下降。

本文計算鈾氧化物熔池凝固殼層的厚度如圖6(a)所示，而文獻[2]根據(jù)嚴重事故下堆芯熔融物在壓力容器內(nèi)滯留的分析程序IVRASA計算得出的數(shù)據(jù)，以及USCB (DOE/ID-10460)和INEEL的實驗數(shù)據(jù)如圖6(b)所示?？梢钥吹剑?1) 凝固殼層厚度在壓力容器下封頭底部（傾角為零）最大，約為0.06m，隨著傾角增大，殼層厚度逐漸減小，殼層頂部厚度約為 0.002m；(2) 本文計算的凝固殼層厚度數(shù)值及隨傾角的變化規(guī)律與USCB (DOE/ID-10460)的實驗結(jié)果一致性較好。

圖5 封頭壁面內(nèi)側(cè)的加熱熱流(a) 本文計算，(b) 文獻[2]計算Fig.5 Heat flux at internal wall of downward part. (a) Our computation, (b) Ref.[2]

圖6 鈾氧化物熔池凝固殼層的厚度(a)，文獻[2]凝固殼層的計算與實驗數(shù)據(jù)(b)Fig.6 Thickness of condensed layer of UO2 (a), computing and experimental result of condensed layer in Ref.[2] (b).

本文計算在堆芯熔融物加熱作用下，壓力容器下封頭壁面溫度隨時間的發(fā)展過程如圖7所示，圖7中共給出了自初始時刻至其后15000s范圍內(nèi)的計算結(jié)果?？梢钥吹剑?1) 從初始時刻開始，壁面溫度逐漸上升，在1000s時金屬熔池周圍的壁面溫度達到1500K，文獻[2]中的圖也表明，燒蝕點附近的壁面溫度為1500K左右；(2) 當金屬熔池周圍的壁面溫度超過壓力容器壁面的熔點溫度，如在2000s時，壁面發(fā)生燒蝕；(3) 熔池的熱量不斷向壁面深處傳遞，壁面溫度全面上升，在3 000s及以后直至9000s的時間內(nèi)，壁面上的燒蝕區(qū)域不斷擴寬并加深；(4) 9000s以后，壓力容器與熔池、外側(cè)安全殼水池之間趨于熱平衡，壁面溫度分布與燒蝕區(qū)域形狀與大小趨于穩(wěn)定，如10000s和15000s的結(jié)果所示，與圖8所示文獻[2]中ABAQUS計算的壁面燒蝕輪廓相符。

圖7 壓力容器下封頭壁面厚度的變化Fig.7 Variation of wall thickness in downward part of pressure vessel.

圖8 文獻[2]中ABAQUS計算的壁面輪廓Fig.8 Wall configuration of wall in Ref.[2].

4 結(jié)語

采用堆芯熔融物的兩層結(jié)構(gòu)模型，計算熔池對壓力容器下封頭壁面的加熱熱流；建立壁面的移動邊界二維平板導熱模型，計算在熔池加熱和安全殼水池冷卻作用下，壓力容器壁面的動態(tài)燒蝕過程；以AP600為研究對象，計算堆芯熔池形成后15000s范圍內(nèi)，壓力容器下封頭壁面溫度及厚度的變化。

研究表明：(1) 鈾氧化物熔池對下封頭壁面的加熱熱流隨傾角而上升，金屬熔融層對壁面的加熱熱流略低；(2) 內(nèi)壁溫度隨時間上升很快，兩熔融層界面處最先出現(xiàn)燒蝕，之后燒蝕深度和寬度增加；(3) 在9000s以后，熔池-壁面-安全殼水池間趨于熱平衡，壁面溫度分布和燒蝕區(qū)基本不再發(fā)生變化。

1 Frepoli C. An overview of westinghouse realistic large break LOCA evaluation model[J]. Science and Technology of Nuclear Installations, 2008, 498737

2 Theofanous T G, Liu C. In-vessel coolability and retention of a core melt[J]. Nuclear Engineering and Design, 1997, 169: 1–48

3 Guan Z H, Yu H X, Jiang G M. Study on cooling model for debris in lower plenum and countermeasures for prevention of focusing effect[J]. Nuclear Engineering and Design, 2008, 29(5): 72–76

4 Fu X L, Yang Y H, Zhou W H, et al. Analysis of process for core melt and debris pool formation in IVR evaluation for CPR1000[J]. Nuclear Engineering and Design, 2010, 31(5): 102–107

5 Luis E H, Mark H A. A diffusion layer model for steam condensation within the AP600 containment[J]. Nuclear Engineering and Design, 1998, 183: 133–150

6 Rafael M J, Kostadin I. Numerical diffusion and the tracking of solute fields in system codes (Part III): application to a boron dilution transient analysis in the AP600[J]. Nuclear Engineering and Design, 1998, 179: 345–356

7 Mark H A, Luis E H. Experimental analysis of heat transfer within the AP600 containment under postulated accident conditions[J]. Nuclear Engineering and Design, 1998, 185: 153–172

8 Li L, Zang X N. The research on molten corium coolability measures for the PWR severe accident[J]. Nuclear Security, 2007, 4: 39–44

9 趙靜, 余紅星, 李鋒. 1000MW級壓水堆安全殼壓力溫度計算分析[J]. 核動力工程, 2003, 24(5): 409–411

ZHAO Jing, YU Hongxing, LI Feng. Calculation and analyses of containment pressure and temperature for PWR 1000MW[J]. Nuclear Engineering and Design, 2003, 24(5): 409–411

10 陶文銓. 數(shù)值傳熱學[M]. 第二版. 西安: 西安交通大學出版社, 2001 TAO Wenquan. Numerical heat transfer[M]. 2ndEd. Xi'an: Xi'an Jiao tong University Press, 2001

11 Globe S, Dropkin D. Natural-convection heat transfer in liquids confined by two horizontal plates and heated from below[J]. Journal of Heat Transfer, 1959, 81: 24–28

CLC TL99

Numerical simulation of the ablation process of the nuclear pressure vessel heated by core melt

ZHANG Xiaoying1YAO Tingting1LI Zhiwei1HUANG Kai2
1(School of Electricity, South China University of Technology, Guangzhou 510640, China) 2(China Nuclear Power Technology Research Institute, Shenzhen 518026, China)

Background: The study on the transient ablation of the pressure vessel wall, acted by the detained melt core, is very important for damage prediction after occurrence of serious loss of coolant accident (LOCA), and design of accident mitigation scheme. Purpose: This paper attempts to study the ablation process of downward part of pressure vessel of AP600, in condition of losing all coolant. Methods: A two-layer of melt material model was used to simulate the structure of melt core, and the heat flux on the vessel wall was computed. The unsteady 2-D heat conduction model was developed for the downward part of pressure vessel, and water boiling in containment was considered as cooling of outer wall. A moving boundary model was used to simulate the decreasing wall thickness on some point in the ablation process. The variation of wall temperature and thickness with time, from occurrence of the accident, till 3000s after the accident, was computed in detail. Results and Conclusion: Simulation results shows that heat flux from the molten UO2on the vessel wall grows up with increasing inclination angle, and heat flux from the molten metal is a little small. Temperature on internal wall grows up rapidly with time, which causes the ablation firstly appearing at interface of the two molten layers, and the ablation zone tends to deepen and widen after that. The heat balance between molten material, vessel wall and containment pool will come 9000s after occurrence of accident, and temperature and ablation zone would not change after 9000s.

Reactor, Coolant loss, Core melt, Pressure vessel, Wall ablation

TL99

10.11889/j.0253-3219.2015.hjs.38.020606

項目(No.51176052、No.51376065)資助

張小英，女，1973年出生，2002年于北京航空航天大學獲博士學位，研究領(lǐng)域為反應堆熱工水力計算、核電安全理論與技術(shù)

黃凱，E-mail: hkai@cgnpc.com.cn

2014-01-12，

2014-05-26