嚴(yán)重事故IVR下反應(yīng)堆壓力容器耦合傳熱數(shù)值模擬分析

姚彥貴，施 楊，蔣 興，賀寅彪

(上海核工程研究設(shè)計(jì)院，上海 200233)

先進(jìn)輕水堆(即第三代)核電廠AP1000設(shè)計(jì)遵循“ALWR-URD”，提高了核電廠安全指標(biāo)的要求，將堆芯熔化物堆內(nèi)包容(IVR)作為重要的緩解嚴(yán)重事故的安全措施之一。IVR是當(dāng)發(fā)生堆芯融化初期，冷卻水通過(guò)RPV外的金屬保溫層底部的注入RPV和金屬保溫層之間通道內(nèi)，冷卻RPV外壁，通過(guò)沸騰換熱帶走RPV下封頭內(nèi)熔融物產(chǎn)生的衰變熱，從而使熔融物包容在RPV內(nèi)，起到緩解事故的作用。因此堆腔注水將RPV傳遞的熱量充分帶走，保持在水沸騰換熱的臨界熱負(fù)荷以內(nèi)，是保證RPV壓力邊界結(jié)構(gòu)完整性的必要條件。而防止RPV隨時(shí)間的蠕變失效則是保證RPV壓力邊界完整性的充要條件。RPV主要存在如下兩個(gè)問(wèn)題：1) 堆腔注水是否保證RPV的有效冷卻與堆內(nèi)熔融物熱量釋放；2) 是否能保證防止RPV的高溫蠕變失效[1-5]。而開(kāi)展嚴(yán)重事故IVR下RPV結(jié)構(gòu)完整性分析的主要輸入是RPV內(nèi)外壁的溫度分布和RPV壁發(fā)生熔融后的剩余壁厚。

近年來(lái)，Sehgal等[6-7]、美國(guó)SNL的Pilch等[8]、法國(guó)CEA的Devos等[9]對(duì)RPV進(jìn)行了一系列的縮比實(shí)驗(yàn)研究。但縮比實(shí)驗(yàn)并不能完全模擬RPV內(nèi)外的自然對(duì)流情況，且實(shí)驗(yàn)代價(jià)昂貴。德國(guó)Willschutz采用有限元方法對(duì)RPV進(jìn)行溫度場(chǎng)和應(yīng)力場(chǎng)分析，并與FOREVER實(shí)驗(yàn)數(shù)據(jù)對(duì)比；但該計(jì)算基于給定RPV壁面邊界條件進(jìn)行，無(wú)法模擬RPV內(nèi)外流體流動(dòng)換熱情況。本文采用流-固-熱耦合的方式，對(duì)RPV內(nèi)熔融物、RPV壁面及RPV外部氣液兩相流進(jìn)行耦合數(shù)值分析，以獲得耦合情況下的溫度場(chǎng)、流場(chǎng)及RPV熔化情況，為RPV在IVR下的結(jié)構(gòu)完整性分析提供重要依據(jù)。

1 數(shù)學(xué)物理模型

IVR-RPV耦合現(xiàn)象屬于氣液兩相流-固體-熔融物耦合換熱流動(dòng)，在該過(guò)程中包含自然對(duì)流沸騰換熱及固體熔化、局部熔融物冷凝等問(wèn)題，為使計(jì)算可行，需進(jìn)行相應(yīng)的簡(jiǎn)化假設(shè)：1) RPV外部堆腔冷卻水為核態(tài)沸騰；2) RPV內(nèi)部熔融物為兩層結(jié)構(gòu)；3) 熔融物底層為氧化層，上層為金屬層；4) 熔融物上部空間空氣為理想氣體；5) RPV內(nèi)壁的熔融過(guò)程通過(guò)迭代計(jì)算完成，即第一步假設(shè)RPV壁厚完整情況下計(jì)算RPV的內(nèi)壁溫度分布，通過(guò)計(jì)算獲得的內(nèi)壁溫度分布，判斷如果大于RPV熔點(diǎn)溫度，表示該部分已熔化，這樣在第二步中重新建模，并將熔化的部分挖掉，這樣反復(fù)迭代幾次，直到RPV內(nèi)壁溫度低于熔點(diǎn)溫度。

本文采用CFD通用商業(yè)軟件STAR-CCM+進(jìn)行耦合傳熱數(shù)值模擬計(jì)算，計(jì)算區(qū)域模型如圖1所示，分析計(jì)算所需數(shù)學(xué)物理模型如下。

圖1 模擬計(jì)算區(qū)域示意圖

1.1 流體域控制方程

計(jì)算區(qū)域的流體域包括RPV外部?jī)上嗔魉?、熔融物、RPV熔化層及熔融物頂部空氣區(qū)域。流體域控制方程包含質(zhì)量、動(dòng)量、能量守恒方程：

ρα·(rαραUα

(1)

(2)

(3)

1.2 固體域控制方程

固體域?qū)岱匠倘缦拢?/p>

λ

(4)

其中：cV為比定容熱容；T為溫度；t為時(shí)間；λ為導(dǎo)熱系數(shù)。

1.3 相間動(dòng)量交換模型

在氣液兩相流動(dòng)過(guò)程中，任一α相會(huì)受到其他相作用力，包括相間拽力和相間非拽力，可表示為：

(5)

(6)

式中：db為氣泡平均直徑，m；CD取決于氣泡雷諾數(shù)Reb，本文采用Ishii-Zuber[10]關(guān)系式計(jì)算；ρf為流體密度；Ug和Uf為氣相和液相速度。

(7)

式中：FTD為湍流耗散力；CTD為湍流耗散力系數(shù)；Ccd為相間拽力動(dòng)量傳遞系數(shù)，由Burns等[11]的經(jīng)驗(yàn)關(guān)系式計(jì)算；vtc為修正系數(shù)；σtc為時(shí)間常數(shù)。

1.4 相間熱量交換模型

相間傳熱采用雙熱阻模型，液相與界面間的傳熱為：

qls=hls(Ts-Tl)

(8)

其中：Ts為交界面溫度；Tl為液相溫度；hls為界面?zhèn)鳠嵯禂?shù)，采用Ranz-Marshall[12]經(jīng)驗(yàn)關(guān)系式計(jì)算。

0≤Reb≤200

(9)

在氣相側(cè)，采用零熱阻模型，其hα為無(wú)窮大，等效于交界面溫度，即氣相飽和溫度。

1.5 相間質(zhì)量交換模型

采用Rohsenow模型[13]計(jì)算蒸發(fā)率，沸騰壁面熱流密度為：

(10)

其中：μl為液相動(dòng)力黏度；hlat為汽化潛熱；ρl、ρv分別為液相、氣相密度；g為重力加速度；σ為表面張力；cpl為液相比定壓熱容；Tw、Tsat分別為壁面和飽和溫度；Cqw為經(jīng)驗(yàn)系數(shù)；Prl為液相普朗特?cái)?shù)。

壁面核化點(diǎn)的蒸發(fā)率為：

(11)

其中，Cew為經(jīng)驗(yàn)系數(shù)，表示壁面熱流用于蒸發(fā)的份額。

2 網(wǎng)格劃分及邊界條件

2.1 網(wǎng)格劃分

考慮到計(jì)算區(qū)域幾何對(duì)稱性，分析采用二維軸對(duì)稱模型，為了模擬冷卻水的自然循環(huán)，模型中建立水槽，并形成一定的水位，水位形成的壓頭為冷卻水的入口壓力，如圖2所示?？疾鞂?duì)象為RPV的熔化，因此對(duì)近壁面區(qū)域進(jìn)行網(wǎng)格加密，對(duì)外水槽區(qū)域加粗。為了進(jìn)行獨(dú)立解驗(yàn)證，共采用3套網(wǎng)格劃分方案，經(jīng)獨(dú)立解分析，確定最終選用的網(wǎng)格數(shù)量為64 023，網(wǎng)格示意圖如圖3所示。

圖2 分析模型

圖3 分析模型網(wǎng)格示意圖

2.2 邊界條件

為確定偏微分方程的解，需設(shè)定合適的邊界條件，本文依據(jù)物理背景，給定邊界條件如下：

1) 氣液兩相流與外壁面邊界：絕熱、無(wú)滑移壁面；

2) 其他液體與固體交界面：熱耦合無(wú)滑移壁面，熱流密度由程序自動(dòng)計(jì)算；

3) 氧化層：給定體積熱源；

4) 水隙頂部邊界：自由液面；

5) 金屬層上表面：輻射換熱；

6) 空氣頂部邊界：壓力入口邊界。

2.3 物理模型

計(jì)算所采用的物理模型有：1) 二維軸對(duì)稱計(jì)算；2) 隱式非定常；3) 多相流模型，包括歐拉多相流模型、VOF方法和沸騰模型；4) 分離式求解；5) 重力模型；6) 湍流模型，即Realizable K-Epsilon模型。

計(jì)算需考慮VOF模型的沸騰模型，因此采用多相流的自由界面模型(VOF方法)，并考慮重力影響。同時(shí)本次VOF計(jì)算需采用非穩(wěn)態(tài)計(jì)算，最終計(jì)算至液膜達(dá)到穩(wěn)定。

2.4 介質(zhì)物性

表1列出沸騰相關(guān)物性參數(shù)，表2列出凝固熔化相關(guān)物性參數(shù)。

表1 沸騰相關(guān)物性參數(shù)

表2 凝固熔化相關(guān)物性參數(shù)

3 模擬結(jié)果與分析

圖4 溫度分布

圖4為溫度分布云圖，可見(jiàn)氧化物區(qū)域溫度較高，且靠近金屬層處溫度最高，達(dá)到3 500 K左右，靠近底部溫度逐漸降低，溫度分布出現(xiàn)明顯的熱分層現(xiàn)象。右圖為根據(jù)溫度分布預(yù)測(cè)的RPV壁熔解區(qū)域，RPV壁材料為SA508 G3鋼，其熔點(diǎn)為1 600 K，根據(jù)溫度分布只要高于1 600 K就會(huì)熔解。從圖中可知，靠近金屬層處還會(huì)出現(xiàn)熔解現(xiàn)象，與金屬層接觸的只有靠近氧化物層的底部出現(xiàn)熔融現(xiàn)象。

圖5 固體體積分率分布

圖5為熔融物區(qū)域固體體積份額分布，可看出，除底部出現(xiàn)大量凝固外，在RPV交界面上也有一薄層凝固層。這是因?yàn)樵赗PV內(nèi)壁面上，雖然此時(shí)RPV被加熱甚至熔化，但是該處溫度仍遠(yuǎn)低于熔融物熔化溫度(凝固溫度)，由于熱交換的作用，靠近RPV內(nèi)壁面的熔融物溫度降低到凝固點(diǎn)之下，從而產(chǎn)生一層凝固層。

圖6為速度矢量分布，可知空氣域存在由中心上升、壁面附近下降的自然對(duì)流，而沸騰域存在周期性運(yùn)動(dòng)，1個(gè)周期內(nèi)存在向上或向下的運(yùn)動(dòng)。熔池內(nèi)熔融物也發(fā)生自然對(duì)流，熔融物沿RPV內(nèi)壁面向下流動(dòng)。圖7為底部水槽進(jìn)口流量隨時(shí)間的變化，可知進(jìn)入沸騰區(qū)域的水發(fā)生周期性的變化，這主要是由于流道內(nèi)兩相流阻力的不穩(wěn)定，與恒定的驅(qū)動(dòng)壓頭之間的不斷平衡造成的。

圖6 速度矢量分布

圖7 進(jìn)口水流量隨時(shí)間的變化

圖8為RPV內(nèi)外壁面溫度及熱流密度分布，x坐標(biāo)為中心軸向坐標(biāo)，中心軸為從RPV底部到金屬層上表面，模型中從RPV到金屬層上表面共2.5 m。從圖8可知，在氧化物層與RPV壁交界的中心高度熱流密度出現(xiàn)峰值，且一直到金屬層交界處，壁面熱流密度達(dá)2 200 kW/m2，這也是金屬層導(dǎo)熱系數(shù)高的原因，因此氧化物層傳給金屬層的熱迅速向RPV壁方向傳遞。從RPV內(nèi)壁面溫度分布可知，內(nèi)壁面溫度高于熔點(diǎn)，因此與氧化物層交界處發(fā)生熔融現(xiàn)象。

a,c:——RPV內(nèi)壁氧化物區(qū)，——RPV內(nèi)壁金屬層區(qū)；b,d:——RPV外壁氧化物區(qū)，——RPV外壁金屬層區(qū)

4 結(jié)論

本文采用大型CFD分析軟件STAR-CCM+，對(duì)嚴(yán)重事故IVR下RPV內(nèi)熔融物和堆外冷卻水沸騰換熱進(jìn)行了耦合數(shù)值模擬分析，通過(guò)耦合分析模擬了嚴(yán)重事故下RPV的溫度分布。為嚴(yán)重事故IVR下RPV結(jié)構(gòu)完整性分析提供了依據(jù)，并得到以下結(jié)論：

1) IVR情況下，熔融物內(nèi)部會(huì)形成熱分層，在熔融物底部與RPV內(nèi)壁面交界處發(fā)生凝固現(xiàn)象；

2) RPV外部的冷卻水會(huì)形成波動(dòng)性的氣液兩相流；

3) RPV內(nèi)壁面部分材料被熔化，RPV壁面未被熔穿，仍保留一定的壁厚；

4) 在極限熱源條件下，RPV熔融后壁厚最薄處達(dá)10 mm，壁面熱流密度最高達(dá)2 200 kW/m2。

本文工作開(kāi)展過(guò)程中得到了艾迪捷信息科技(上海)有限公司軟件技術(shù)人員的大力支持，在此對(duì)李晶博士、江興賢和韓海博士表示衷心的感謝。

參考文獻(xiàn)：

[1] SEHGAL B R, THEERTHAN A, GIRI A, et al. Assessment of reactor vessel integrity (ARVI)[J]. Nuclear Engineering and Design, 2003, 221(1-3): 23-53.

[2] 姚偉達(dá)，謝永誠(chéng)，賀寅彪，等. 第三代非能動(dòng)型與第四代超臨界先進(jìn)輕水型核電廠中典型的力學(xué)問(wèn)題[C]∥第14屆全國(guó)反應(yīng)堆結(jié)構(gòu)力學(xué)會(huì)議論文集. 桂林：中國(guó)力學(xué)學(xué)會(huì)，2006.

[3] ASME boiler and pressure vessel code: Rules for construction of nuclear facility components[S]. New York, USA: ASME, 2004.

[4] 武志瑋，寧冬，姚偉達(dá). 嚴(yán)重事故下反應(yīng)堆壓力容器材料高溫蠕變研究進(jìn)展[J]. 核安全，2011(2)：20-24.

WU Zhiwei, NING Dong, YAO Weida. Research progress on high-temperature creep behavior of reactor pressure vessel[J]. Nuclear Safety, 2011(2): 20-24(in Chinese).

[5] 姚彥貴，寧冬，武志瑋，等. 假想堆芯熔化嚴(yán)重事故下反應(yīng)堆壓力容器完整性的研究進(jìn)展與建議[J]. 核技術(shù)，2013，36(4)：040615-1-040615-5.

YAO Yangui, NING Dong, WU Zhiwei, et al. Research progress and recommendations on reactor pressure vessel integrity under hypothetical core melt down accident[J]. Nuclear Techniques, 2013, 36(4): 040615-1-040615-5(in Chinese).

[6] SEHGAL B R, NOURGALIEV R R, DINH T N. Characterization of heat transfer processes in a melt pool convection and vessel-creep experiment[J]. Nuclear Engineering and Design, 2002, 211(2-3): 173-187.

[7] WILLSCHüTZ H G, ALTSTADT E, SEHGAL B R, et al. Coupled thermal structural analysis of LWR vessel creep failure experiments[J]. Nuclear Engineering and Design, 2001, 208(3): 265-282.

[8] PILCH M M, LUDWIGSEN J S, CHU T Y, et al. Creep failure of a reactor pressure vessel lower head under severe accident conditions[C]∥ASME/JSME Joint Pressure Vessel and Piping (PVP) Conference. San Diego, CA: ASME, 1998.

[9] DEVOS J, SAINTE C C, POETTE C, et al. CEA programme to model the failure of the lower head in severe accidents[J]. Nuclear Engineering and Design, 1999, 191(1): 3-15.

[10] ISHII M, ZUBER N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows[J]. AIChE J, 1979, 25: 843-855.

[11] BURNS A D, FRANK T, HAMILL I, et al. The FAVRE averaged drag model for turbulent dispersion in Eulerian multi-phase flows[C]∥5th International Conference on Multiphase Flow, ICMF-2004. Yokohama, Japan: JSME, 2004.

[12] RANZ W E, MARSHALL W R. Evaporation form drops[J]. Chem Eng Prog, 1952, 48(3): 141-147.

[13] ROHSENOW W M. A method of correlation heat transfer data for surface boiling of liquid[J]. Trans ASME, 1952, 74: 969-978.