反應(yīng)堆內(nèi)熔融物冷卻的三維數(shù)值模擬研究

薛 峰，袁明豪，張 建，陳秋煬

(蘇州熱工研究院有限公司 核安全與運行技術(shù)中心，江蘇 蘇州 215004)

核電廠發(fā)生嚴(yán)重事故，特別是在發(fā)生高壓熔堆事故時，會對安全殼的完整性造成極大威脅，因此如何保持壓力容器完整性，將堆芯熔融物保持在壓力容器內(nèi)是嚴(yán)重事故緩解措施研究的重點。堆芯熔融物壓力容器內(nèi)滯留(IVR)是一種重要的嚴(yán)重事故管理措施，主要應(yīng)用于先進(jìn)的輕水堆核電站，其主要目的為通過壓力容器外部的冷卻降低內(nèi)部熔融物的溫度，從而保持壓力容器的完整性，使熔融物滯留在壓力容器內(nèi)部，避免放射性物質(zhì)的釋放。

國內(nèi)外對IVR進(jìn)行了大量研究和應(yīng)用。芬蘭Lovvisa核電廠首次在VVER-440反應(yīng)堆上應(yīng)用了堆腔注水系統(tǒng)[1]。西屋的AP600和AP1000[2]以及韓國的APR1400[3]在設(shè)計中也考慮了堆腔注水概念，開展了相關(guān)的實驗和計算驗證其安全性。中國自主設(shè)計的華龍一號中已應(yīng)用IVR技術(shù)[4]，以提高核電廠應(yīng)對嚴(yán)重事故的能力。

堆芯發(fā)生熔化后，熔融物坍塌至壓力容器下封頭形成熔池，熔池自身的特性和結(jié)構(gòu)對于壁面?zhèn)鳠岷拖路忸^的完整性具有重要的影響。楊曉等[5]開發(fā)了MOPOL程序并對IVR有效性評價中的不確定性進(jìn)行了分析。陳星等[6]和傅孝良[7]分別使用MOPOL程序?qū)PR1000堆型的IVR有效性進(jìn)行了評價。Zhang等[8-9]開發(fā)了IVR分析的軟件IVRASA并對AP1000的不同工況進(jìn)行了評價分析，另外還通過實驗裝置COPRA開展了下封頭內(nèi)堆芯熔融物換熱特性的實驗研究。曹臻等[10]開發(fā)了三層熔池模型分析程序SPIRE并進(jìn)行了驗證。本文基于日本應(yīng)用能源研究所(IAE)開發(fā)的SAMPSON程序，對壓力容器內(nèi)熔融物擴(kuò)展和冷卻分析(DCA)模塊進(jìn)行改進(jìn)開發(fā)，并與實驗結(jié)果進(jìn)行對比。

1 模型及基本方程

1.1 熔融物擴(kuò)展冷卻模型

SAMPSON程序DCA模塊的瞬態(tài)計算流程如圖1所示。

圖1 DCA模塊的瞬態(tài)計算流程

熔融物在壓力容器下腔室內(nèi)的擴(kuò)展和冷卻模型如圖2所示，該模型的分析方法如圖3所示。連續(xù)相的融熔物被視為牛頓流體，流動和傳熱采用N-S方程進(jìn)行分析。網(wǎng)格被賦予兩種屬性：自然對流分析網(wǎng)格及液相自由表面網(wǎng)格，以模擬熔融物的擴(kuò)展、熔化和凝固的作用。自然對流分析網(wǎng)格是指充滿融熔物、硬殼的網(wǎng)格，采用SMAC(simplified marker and cell)方法進(jìn)行流動計算。液相自由表面網(wǎng)格是指沒有被充滿的網(wǎng)格，采用高度函數(shù)來計算流動。

圖2 熔融物在壓力容器下腔室內(nèi)的擴(kuò)展和冷卻模型

本文分析中作如下假設(shè)：1) 連續(xù)相融熔物的擴(kuò)展是由重力主導(dǎo)的，因此不考慮融熔物向空氣散布及其中涉及的氣相；2) 硬殼的遷移被限制在豎直方向。由假設(shè)1可知，從自然對流分析網(wǎng)格流到分析邊界外面網(wǎng)格的融熔物會聚集在液相自由表面網(wǎng)格的下部。

圖3 熔融物擴(kuò)展和冷卻模型的分析方法

1.2 基本方程

1) 自然對流分析網(wǎng)格的基本方程

DCA模塊中關(guān)于自然對流分析網(wǎng)格內(nèi)的擴(kuò)展和自然對流行為的基本方程包括質(zhì)量守恒、動量守恒和能量守恒方程。

質(zhì)量守恒方程為：

(1)

動量守恒方程為：

(2)

(3)

(4)

其中：u、v、w是流體在t時刻，在坐標(biāo)點(x、y、z)處的速度分量；ρ為密度；p為壓力；μ為動力黏度；g為重力加速度；H為高度；K為浮升力項，采用Boussinesq假設(shè)。

(5)

能量守恒方程為：

(6)

其中：e為比焓；λ為導(dǎo)熱系數(shù)；T為溫度；Q為熱源。

熔融物的熔化、凝固狀態(tài)，由以下公式計算每個網(wǎng)格的溫度和凝固份額來確定，其中比焓由能量守恒方程得到。

(7)

其中：cs和cp分別為固態(tài)和液態(tài)時的比熱容；Tl為熔化溫度；el為熔化比焓；es為凝固比焓；b為凝固份額，完全凝固時，b=1，完全熔化時，b=0。當(dāng)b超過某一值時，流動將停止。對于鐵和非鐵合金，b在0.5～0.7范圍內(nèi)。

當(dāng)b超過流動限值時，假設(shè)連續(xù)相融熔物凝固，流體速度被置為0。通過在每個計算步長的壓力計算中應(yīng)用該邊界條件，來實現(xiàn)流動計算中對固液界面遷移的處理。對于能量計算，由于可計算凝固網(wǎng)格的導(dǎo)熱，固相和液相被當(dāng)作連續(xù)體。

2) 液相自由表面網(wǎng)格的基本方程

根據(jù)假設(shè)，連續(xù)相融熔物的遷移受到重力和壁面的限制。對于液位參考平面不變的情況，使用高度函數(shù)來追蹤液相表面的位置。xy平面作為參考平面，z方向表示距離參考平面的液位。液體高度H的運動方程為：

(8)

1.3 數(shù)值算法

DCA模塊的流場計算采用SMAC方法。首先求解自然對流分析邊界內(nèi)的質(zhì)量守恒和動量守恒方程，得到連續(xù)相融熔物的流速，然后使用該流體速度求解能量方程。

網(wǎng)格采用固定交錯網(wǎng)格，壓力和溫度定義在網(wǎng)格的中心，流速定義在網(wǎng)格的界面上，如圖4所示。圖4中，i、j、k分別為x、y、z方向上的網(wǎng)格編號。

2 程序開發(fā)

2.1 兩層熔池模型

熔池在壓力容器下封頭達(dá)到最終穩(wěn)定分層結(jié)構(gòu)之前要經(jīng)歷的中間過程仍存在較大的不確定性，國內(nèi)外[7,11]提出了不同的熔池分層結(jié)構(gòu)，本文采用目前最具有代表性的兩層熔池模型。熔池底層為密度較大的氧化池，包含UO2、ZrO2等氧化物和內(nèi)熱源。此外，氧化池與溫度較低的下封頭和上部金屬層接觸，冷凝成一層硬殼。由于硬殼頂部覆蓋的金屬層溫度與下封頭容器壁溫度不同，故氧化池上、下硬殼厚度并不一致。

a——壓力、焓計算二維網(wǎng)格；b——速度計算二維網(wǎng)格圖4 交錯網(wǎng)格示意圖

熔池上層由密度相對較小的金屬層構(gòu)成，包含熔化的不銹鋼和未被氧化的鋯等金屬物。由于金屬層質(zhì)量與氧化池質(zhì)量相比往往很小，故其厚度較薄。

無論采用何種熔池結(jié)構(gòu)，當(dāng)堆芯熔化并塌落至下封頭形成分層熔池后，下封頭壁面受熱的熱流密度沿其軸向是非均勻分布的。熔池上層的金屬層由于質(zhì)量小、厚度薄，造成其熱流密度很大，形成熱聚集效應(yīng)。若此處的熱流密度超過臨界熱流密度(CHF)，會對壓力容器的完整性構(gòu)成巨大威脅。

在SAMPSON程序的DCA模塊中，假設(shè)熔融物的成分是均勻的，無法模擬熔融物中的金屬層的聚焦作用。為此本文對該模塊的下封頭內(nèi)熔融物冷卻模型進(jìn)行改進(jìn)，增加氧化池與金屬分層模型，如圖5所示。

圖5 氧化池與金屬分層模型

在熔融物擴(kuò)展和自然對流模擬過程中，假設(shè)金屬與氧化池分層布置，各自的厚度(體積)取決于掉落到下封頭內(nèi)的熔融物的質(zhì)量和組分。

金屬層與氧化池的物性分別計算，采用物性包MATPRO[12]。其中金屬的主要成分為鐵和鋯，鋯金屬的密度使用CDEN函數(shù)計算：

ρ=ρ0exp(-εx)exp(-εy)exp(-εz)

(9)

其中：ρ0為參考密度，取溫度為300 K時的值6.55×103kg/m3；εx、εy、εz為任意正交坐標(biāo)系下的應(yīng)變，由CTHEXP函數(shù)計算。因為熱應(yīng)變小于1，因此式(9)可近似為：

ρ=ρ0(1-εx-εy-εz)

(10)

鋯金屬的導(dǎo)熱系數(shù)λz使用CTHCON函數(shù)計算，僅考慮溫度為決定因素，當(dāng)溫度低于2 098 K時，其計算公式為：

λz=7.51+2.09×10-2T+

1.45×10-5T2+7.67×10-9T3

(11)

鋯金屬導(dǎo)熱系數(shù)的實驗數(shù)據(jù)及擬合曲線如圖6所示。

鋼的密度使用SDEN函數(shù)計算，參考密度取溫度為300 K時的值7.9×103kg/m3。

鋼的導(dǎo)熱系數(shù)使用STHCON函數(shù)計算，其公式為：

λs=7.58+0.189T300 K≤T<1 671 K

λs=610.93+0.342T

1 671 K≤T<1 727 K

λs=20 W/(cm·K)T≥1 727 K

(12)

在本文計算中僅考慮鋼完全熔化的工況，因此取其導(dǎo)熱系數(shù)為20 W/(cm·K)。

圖6 鋯金屬的導(dǎo)熱系數(shù)

金屬層與氧化池分界面處的物性采用插值計算，圖7給出了分層界面附近的導(dǎo)熱系數(shù)插值示意圖。

圖7 分層界面附近的導(dǎo)熱系數(shù)插值

假設(shè)氧化池的導(dǎo)熱系數(shù)為λo，金屬層的導(dǎo)熱系數(shù)為λm，界面附近的導(dǎo)熱系數(shù)可采用下式計算：

(13)

λu=λof+λm(1-f)

(14)

其中：Δz為網(wǎng)格中心距；Δz-、Δz+、Δz1、Δz2的含義如圖7所示；λd和λu為包含分界面的網(wǎng)格單元的導(dǎo)熱系數(shù)，分別對應(yīng)分界面位于網(wǎng)格中心的下方和上方(圖7)兩種情況；f為插值系數(shù)。

2.2 壓力容器外表面?zhèn)鳠?/h3>

壓力容器的傳熱模型和網(wǎng)格劃分如圖8所示，其中，qin為來自熔融物的熱流，qcond為導(dǎo)熱熱流，qout為到安全殼大氣的熱流。通過求解三維導(dǎo)熱方程(式(15))得到壓力容器的溫度分布。

圖8 壓力容器壁面網(wǎng)格和傳熱模型

e=ct+(1-b)Hsl

(15)

其中：λ為導(dǎo)熱系數(shù)；Q為傳入熱量；c為比熱容；Hsl為熔化潛熱。

壓力容器下封頭外表面向堆腔的傳熱包含向空氣和水兩部分的傳熱，其計算公式如下：

qrpv=hgas(1-α)Arpv(T(θ)-Tgas)+

hpoolαArpv(T(θ)-Tpool)

(16)

其中：qrpv為下封頭外表面向堆腔的傳熱；hgas和hpool分別為下封頭到堆腔空氣和水的對流傳熱系數(shù)；α為下封頭被水淹沒的面積份額；Arpv為下封頭表面積；T(θ)為下封頭外表面在角度θ處的溫度，其中θ=0°表示下封頭底部；Tgas和Tpool分別為堆腔內(nèi)空氣和水的溫度。

對于外表面與空氣的傳熱，hgas使用大空間自然對流換熱實驗關(guān)聯(lián)式(式(17))計算，此關(guān)聯(lián)式在工程中廣泛使用。

Nu=C(GrPr)n

(17)

其中：Gr為格拉曉夫準(zhǔn)則數(shù)；Nu為努塞爾數(shù)；Pr為普朗特數(shù)；h為對流換熱系數(shù)；L為傳熱面的幾何特征長度；系數(shù)C和指數(shù)n根據(jù)流體處于層流或湍流取不同的值。

外表面與水的傳熱在下封頭外表面溫度未達(dá)到堆腔中水的飽和溫度時，使用式(16)計算自然對流換熱。當(dāng)下封頭外表面溫度大于飽和溫度時，則考慮沸騰傳熱。

2.3 網(wǎng)格前處理程序

由于下封頭內(nèi)的熔融物采用三維直角坐標(biāo)網(wǎng)格，而下封頭本身采用三維曲線坐標(biāo)網(wǎng)格(圖9)，所以必須預(yù)先確定兩套網(wǎng)格之間交界面的幾何參數(shù)。SAMPSON程序的DCA模塊本身不具備產(chǎn)生這些幾何參數(shù)的能力，這些參數(shù)需人工計算產(chǎn)生，當(dāng)網(wǎng)格數(shù)量較大時，產(chǎn)生這些參數(shù)的工作量很大，無法確保其可靠性。因此，為便于快速準(zhǔn)確地得到兩套網(wǎng)格間的幾何參數(shù)，以及互相映射的關(guān)系，本文開發(fā)了前處理程序來自動產(chǎn)生這些參數(shù)。

圖9 壓力容器網(wǎng)格(紅色)和熔融物網(wǎng)格(灰色)

該前處理程序首先根據(jù)設(shè)定的壓力容器下封頭節(jié)點布置產(chǎn)生網(wǎng)格，再將下封頭網(wǎng)格的任一內(nèi)表面分為M×N個極小的微面積，根據(jù)每一微面積的坐標(biāo)判斷其所處的熔融物三維直角坐標(biāo)網(wǎng)格單元。通過對M×N個微面積進(jìn)行累加，即可確定該下封頭網(wǎng)格內(nèi)表面所映射的熔融物三維直角坐標(biāo)網(wǎng)格(可對應(yīng)多個)以及相互接觸的面積。對下封頭內(nèi)表面所有網(wǎng)格循環(huán)計算后，就可得到關(guān)于熔融物網(wǎng)格與壓力容器網(wǎng)格交界面的所有幾何參數(shù)。

3 驗證計算

3.1 計算模型與假設(shè)

為對改進(jìn)后的DCA模塊進(jìn)行驗證，本文選取AP1000核電機(jī)組建立模型進(jìn)行計算分析。AP1000核電機(jī)組下封頭的尺寸根據(jù)文獻(xiàn)[11]選取。壓力容器下封頭網(wǎng)格與熔融物網(wǎng)格的敏感性分析如圖10所示，最終壓力容器下封頭選取42×72×20個網(wǎng)格，壁面方向的網(wǎng)格加密是為了盡可能詳細(xì)模擬下封頭壁面的減薄效應(yīng)。下封頭內(nèi)部的熔融物網(wǎng)格為42×42×32的直角坐標(biāo)網(wǎng)格。熔融物網(wǎng)格與壓力容器網(wǎng)格交界面的映射關(guān)系由前處理程序自動產(chǎn)生。

文獻(xiàn)[11]采用拉丁超立方體抽樣(LHS)技術(shù)對堆芯熔融物材料分布和余熱進(jìn)行計算，給出了基準(zhǔn)工況假設(shè)，本文采用此工況進(jìn)行計算，基準(zhǔn)工況參數(shù)取值列于表1。

堆芯熔融物逐步塌落至下封頭中形成熔池。先期坍塌的熔融物在下封頭中經(jīng)歷一段時間后形成具有氧化池和金屬層的兩層結(jié)構(gòu)，隨后熔化的熔融物經(jīng)歷一段時間后分解并溶入相應(yīng)各層，最終熔池穩(wěn)定于兩層結(jié)構(gòu)。因此熔融物落下的瞬態(tài)過程對下封頭失效的影響是可忽略的。

圖10 熔融物(a)與壓力容器(b)的網(wǎng)格敏感性分析

表1 基準(zhǔn)工況的參數(shù)

在使用本文程序進(jìn)行計算時，對熔融物落下的瞬態(tài)做了簡單化處理。假設(shè)表1中的熔融物以一定流量快速落入下封頭，當(dāng)這些熔融物全部進(jìn)入下封頭后，氧化物與金屬開始分層。程序的收斂準(zhǔn)則為熔融物產(chǎn)生的衰變余熱與從金屬層上部輻射傳出的熱量及通過下封頭外壁面?zhèn)鞒龅臒崃窟_(dá)到平衡(圖8)，熱平衡公式為：

Q衰變熱=Q金屬層輻射+Q壁面熱導(dǎo)出

(18)

當(dāng)金屬熔池周圍的壁面溫度超過壓力容器壁面的熔點溫度(假設(shè)為1 700 K)時，壓力容器下封頭壁面發(fā)生燒蝕。

3.2 計算結(jié)果

圖11示出中心對稱截面上熔融物的流場和溫度分布。當(dāng)下封頭的熔池達(dá)到穩(wěn)態(tài)后，上層的金屬層溫度約為2 020 K，下層的熔融物溫度達(dá)到3 000 K。氧化池的熱量通過金屬層向上和側(cè)面?zhèn)鬟f。金屬層向上通過輻射傳熱的熱量較小，向容器壁的側(cè)向傳熱通道較為暢通，大部分熱量通過金屬層側(cè)面對流傳導(dǎo)給壓力容器壁面。金屬層的導(dǎo)熱性能遠(yuǎn)高于氧化池，能快速導(dǎo)出熱量，因此其溫度要低于氧化池。

圖11 中心截面上熔融物的流場和溫度分布

由圖11的流場可看出，氧化池同時被壓力容器壁面和上層金屬層冷卻，故其自然對流會形成兩個對稱的渦結(jié)構(gòu)，同樣的壓力容器壁面會較大程度地冷卻金屬層，金屬層也會形成兩個對稱的渦結(jié)構(gòu)。

圖12示出壓力容器下封頭壁面和厚度方向的溫度分布。由于金屬層會向壓力容器壁面?zhèn)鬟f較多的熱量，與金屬層接觸的壁面有較高的溫度且會達(dá)到其假設(shè)熔點1 700 K，部分厚度的壁面會被熔化，而氧化池接觸的壓力容器壁面溫度分布較為均勻。故可看出金屬層的特性會對壓力容器的失效起著決定作用。

圖13示出壓力容器下封頭外壁面的熱流密度和被熔化的厚度。由圖13可知，在金屬層與下封頭接觸之處，下封頭厚度急劇減薄，此處的熱流密度急劇升高，金屬層的聚焦作用較為明顯。將圖13的結(jié)果與文獻(xiàn)[11]的進(jìn)行比較，在金屬層與下封頭接觸處的熱流密度兩者較接近，但氧化池的熱流密度分布兩者有較明顯區(qū)別。本文計算得到的氧化池的熱流密度隨角度的變化較為平緩，文獻(xiàn)[11]給出的氧化池的熱流密度在角度大于30°時有顯著上升。這是由于文獻(xiàn)[11]在計算時使用了文獻(xiàn)[13]給出的關(guān)系式以及文獻(xiàn)[2]給出的mini-ACOPO實驗數(shù)據(jù)。

圖12 壓力容器下封頭壁面(a)和厚度方向(b)的溫度分布

圖13 壓力容器下封頭外表面的熱流密度(a)和被熔化的厚度(b)

在mini-ACOPO實驗中，熔融物的上表面未考慮金屬層冷卻，而在實際的堆內(nèi)熔融物冷卻滯留分析中，壓力容器下封頭內(nèi)的氧化池的上方覆蓋了一層較低溫度的金屬層，氧化池上方存在較強(qiáng)換熱，與本文計算工況不符，故本文參考BALI實驗[14]結(jié)果進(jìn)行對比分析。

BALI實驗對氧化池上方金屬層的冷卻作用進(jìn)行了研究，實驗結(jié)果表明，存在頂部冷卻時，60°以上區(qū)域的熱流密度變化較均勻。本文重新計算了BALI實驗工況，并對氧化池的熱流密度進(jìn)行無量綱化處理，將之與實驗結(jié)果的擬合曲線進(jìn)行比較，如圖14所示，兩者吻合較好。因此，由于頂部金屬層的冷卻使得靠近金屬層的氧化池的壁面換熱更加均勻，而遠(yuǎn)離金屬層部分的換熱會隨角度有迅速變化的趨勢。

4 結(jié)論

本文對日本應(yīng)用能源研究所(IAE)開發(fā)的SAMPSON程序進(jìn)行二次開發(fā)，在DCA模塊中增加了下封頭內(nèi)堆芯熔融物金屬和氧化物分層模型，完善了壓力容器外表面?zhèn)鳠崮Ｐ?，并獨立開發(fā)了直角坐標(biāo)和曲面坐標(biāo)網(wǎng)格前處理程序。

圖14 下封頭外表面熱流密度與BALI實驗結(jié)果對比

為驗證改進(jìn)的DCA模塊，選取AP1000核電機(jī)組建立下封頭模型進(jìn)行計算分析，并與文獻(xiàn)和實驗結(jié)果進(jìn)行對比，結(jié)果吻合較好，表明改進(jìn)的DCA模塊可有效地模擬下封頭內(nèi)堆芯熔融物的流動和冷卻特性，并用于堆內(nèi)熔融物冷卻滯留的分析評價。