朱光昱 郭 超 劉巧鳳 李 春 依 巖

1（生態(tài)環(huán)境部核與輻射安全中心 北京 100082）

2（中國(guó)核電工程有限公司 核電安全研究中心 北京 100840）

堆芯捕集器是核電廠嚴(yán)重事故緩解手段之一，該技術(shù)通過(guò)在反應(yīng)堆壓力容器（Reactor Pressure Vessel，RPV）外部設(shè)置的收集裝置對(duì)熔融物進(jìn)行冷卻并最終實(shí)現(xiàn)滯留。目前，國(guó)內(nèi)外采用堆芯捕集器設(shè)計(jì)包括：俄羅斯OKB Gidropress 設(shè)計(jì)開(kāi)發(fā)的VVER［1-2］（Vodo-Vodyanoi Energetichesky Reactor）；由法國(guó)和德國(guó)聯(lián)合開(kāi)發(fā)的改進(jìn)型壓水堆EPR［3-4］（Europe Pressurized Water Reactor）；由通用電氣公司和東芝聯(lián)合設(shè)計(jì)的新型沸水堆ESBWR（The Economic Simplified Boiling Water Reactor）［5］；由韓國(guó)電力公司設(shè)計(jì)的APR1400 及其歐洲版本EUAPR1400（Europe-Advanced Power Reactor）［6］。

上述設(shè)計(jì)中，VVER和ESBWR采用了坩堝式設(shè)計(jì)而EPR和APR1400采用了熔融物擴(kuò)展冷卻方式。其中，VVER 的堆芯捕集器［1-2］主體布置在壓力容器下方的堆坑中，坩堝外壁面設(shè)置有冷卻水通道。在事故發(fā)生后，首先通過(guò)犧牲材料與堆芯熔融物發(fā)生反應(yīng)，降低熔融物整體的體積釋熱率并改善熔融池流動(dòng)性。我國(guó)學(xué)者結(jié)合上述設(shè)計(jì)提出了一種雙層坩堝形式的堆芯捕集器，通過(guò)增設(shè)冷卻管提升熔融物冷卻速率［7］。

EPR 堆芯捕集器［3-4］設(shè)計(jì)采用熔融物擴(kuò)展的方式增加換熱面積，從而使得堆芯熔融物快速冷卻。事故發(fā)生后首先通過(guò)堆腔內(nèi)的熔融物臨時(shí)滯留裝置，將分批轉(zhuǎn)移的熔融物進(jìn)行收集，再統(tǒng)一將熔融物鋪展在面積約為170m2的捕集器擴(kuò)展室中，隨后通過(guò)捕集器底面、側(cè)壁面，以及熔融物頂部淹沒(méi)的方式對(duì)熔融物進(jìn)行冷卻。EPR堆芯捕集器預(yù)計(jì)可將熔融物的時(shí)間縮短至幾天之內(nèi)，但系統(tǒng)運(yùn)行流程較多，熔融物擴(kuò)展流動(dòng)的不確定性對(duì)系統(tǒng)設(shè)計(jì)功能的可靠性影響較大，而且堆芯捕集器擴(kuò)展室占據(jù)了大量的空間不利于安全殼內(nèi)的設(shè)備布置。

為了研究坩堝式堆芯捕集器中熔融物冷卻過(guò)程，本文采用多物理場(chǎng)耦合計(jì)算軟件COMSOL，結(jié)合VVER 堆芯捕集器設(shè)計(jì)資料搭建計(jì)算模型，對(duì)熔融物冷卻過(guò)程中的熔融池的溫度、結(jié)殼分布等參數(shù)進(jìn)行了數(shù)值模擬研究，探討了事故緩解進(jìn)程，為該系統(tǒng)的安全性評(píng)價(jià)提供了一定的依據(jù)。

2 計(jì)算模型設(shè)置

2.1 幾何與邊界條件

VVER 堆芯捕集器整體結(jié)構(gòu)如圖1所示，主要由位于RPV 下部的托盤(pán)和位于堆腔內(nèi)的換熱器組成。托盤(pán)內(nèi)設(shè)置有肋片以協(xié)助收集RPV 破口流出的熔融物，并引導(dǎo)熔融物向下部的換熱器遷移。換熱器的內(nèi)部設(shè)置有填料籃來(lái)裝載犧牲材料，其外部設(shè)置有冷卻水通道。

圖1 VVER堆芯捕集器結(jié)構(gòu)示意圖Fig.1 Structure diagram of VVER core catcher

圖2 為參考VVER堆芯捕集器的換熱器部分建立的軸對(duì)稱計(jì)算域模型。坩堝式換熱器主體為圓柱體，半徑約2.355m，底部為圓錐體，錐段與水平方向的傾角為12°。根據(jù)鄭濱［8］對(duì)VVER 堆芯捕集器內(nèi)熔融池的計(jì)算結(jié)果，在分層階段氧化物層和金屬層體積分別為26.1m3和4.6m3，由此得到模型中氧化物層高度為1.61m，金屬層在中軸處的高度為0.5m。熔融物的衰變熱全部集中在氧化物層中，事故發(fā)生后熔融池內(nèi)的衰變熱變化在文獻(xiàn)［8］中有詳細(xì)介紹。

圖2 VVER堆芯捕集器的換熱器軸對(duì)稱計(jì)算域模型示意圖Fig.2 Schematic of axisymmetric computational domain model for heat exchanger of VVER core catcher

計(jì)算過(guò)程中，考慮安全殼內(nèi)壓力維持在0.12MPa，在系統(tǒng)長(zhǎng)時(shí)間運(yùn)行后認(rèn)為冷卻水過(guò)冷度很低，因此將坩堝外壁面溫度設(shè)置為接近飽和溫度的378K。通過(guò)附加薄層的方式考慮不銹鋼坩堝對(duì)導(dǎo)熱的影響，根據(jù)設(shè)計(jì)參數(shù)薄層厚度可設(shè)置為60mm。氧化物層頂部與冷卻水直接接觸可認(rèn)為長(zhǎng)期處于沸騰換熱狀態(tài)，因此頂部設(shè)置為對(duì)流換熱條件。由于缺乏實(shí)際參數(shù)，換熱系數(shù)根據(jù)大容積沸騰米海耶夫公式估算［9］：

式中：h為沸騰換熱系數(shù)，W·(m2·K)-1；Tup和Tsat分別表示頂部邊界網(wǎng)格溫度和飽和溫度，K；P為絕對(duì)壓力，Pa。

2.2 模型設(shè)置

通過(guò)COMSOL 軟件相變模型計(jì)算熔融物的熔化凝固過(guò)程，該模型考慮物質(zhì)相變過(guò)程中固液相線溫度之間為糊狀區(qū)，通過(guò)插值方式計(jì)算糊狀區(qū)的密度、比熱容、導(dǎo)熱系數(shù)。在VVER堆芯捕集器熔融池中，金屬層中98%的質(zhì)量為Fe，其他為少量的金屬Al、U 和Zr［8］，因此，層內(nèi)物質(zhì)的相變潛熱和膨脹系數(shù)均設(shè)置為Fe 的物性。金屬層相變溫度區(qū)間也參照Fe的熔點(diǎn)設(shè)置，同時(shí)考慮了混合物熔點(diǎn)會(huì)適當(dāng)降低。氧化物層中熔融物的固液相線間溫差較大，在凝固過(guò)程中黏度會(huì)逐漸變化，因此對(duì)流動(dòng)過(guò)程影響較大。EPR堆芯捕集器設(shè)計(jì)過(guò)程中對(duì)熔融物擴(kuò)展流動(dòng)特性進(jìn)行了大量實(shí)驗(yàn)研究［10-13］，實(shí)驗(yàn)中各元素質(zhì)量分?jǐn)?shù)和氧化物相的固液相線溫度如表1所示，表中同時(shí)給出了VVER氧化物層各元素質(zhì)量分?jǐn)?shù)。

表1 熔融物成分和固液相線溫度Table 1Melt composition and solidus-liquidus temperature

通過(guò)對(duì)比，本文選擇與VVER 氧化物層主要成分接近的VULCANO 實(shí)驗(yàn)參數(shù)設(shè)置氧化物層的相變潛熱和固液相線溫度。湍流場(chǎng)計(jì)算采用低雷諾數(shù)k-ε模型，固液相線間氧化物黏度μ的變化可根據(jù)懸浮液黏度公式（2）計(jì)算：

式中：μl為熔融物液相黏度，Pa·s；θs為固相體積分?jǐn)?shù)，可直接在計(jì)算過(guò)程中調(diào)用；系數(shù)C取值范圍為3.5～8。

計(jì)算域內(nèi)固相體積分?jǐn)?shù)可以直接從相變模型中實(shí)時(shí)調(diào)取。Journeau等［13］為了簡(jiǎn)化計(jì)算提出了將黏度與溫度擬合的估算方式，同時(shí)推薦系數(shù)C取6，據(jù)此可以將式（2）推導(dǎo)為式（3）：

式中：T、Tliquidus和Tsolidus分別為熔融物溫度、液相線溫度和固相線溫度，K。

文獻(xiàn)［8］中給出了其他計(jì)算所需的參數(shù)，其中采用氧化物層膨脹系數(shù)為0.000 1K-1。在本文中，通過(guò)Boussinesq 假設(shè)計(jì)算熔融池中的自然流動(dòng)過(guò)程，氧化層膨脹系數(shù)參考一般陶瓷類材料的膨脹系數(shù)設(shè)置，由此計(jì)算得到的氧化物層自然循環(huán)速度相對(duì)文獻(xiàn)［8］更加保守。綜合以上內(nèi)容，模擬過(guò)程中使用的物性參數(shù)具體如表2所示。

表2 物性參數(shù)Table 2 Physical properties

2.3 模型參數(shù)敏感性分析

COMSOL Multiphysics 默認(rèn)采用自適應(yīng)時(shí)間步長(zhǎng)來(lái)求解瞬態(tài)問(wèn)題，可以通過(guò)調(diào)整容差因子（Tolerance Factor）可以改變軟件對(duì)于步長(zhǎng)的限制，從而提升計(jì)算效率。模型敏感性分析過(guò)程中，采用事故后2～5 h的衰變熱作為氧化物層加熱功率。

圖3所示為不同網(wǎng)格數(shù)和容差因子下的氧化物層中軸溫度分布，網(wǎng)格數(shù)18 276 和14 472的計(jì)算溫度偏差已經(jīng)縮小至1%以內(nèi)，因此，后續(xù)長(zhǎng)期冷卻模擬采用的網(wǎng)格數(shù)為14 472。通過(guò)在同一臺(tái)計(jì)算機(jī)上對(duì)比研究發(fā)現(xiàn)，網(wǎng)格數(shù)18 276、14 472和11 790的模型中，將容差因子由0.5提高至2可分別將數(shù)值模擬時(shí)間縮短10.86%、10.05%和25.15%，可見(jiàn)，在計(jì)算模型達(dá)到網(wǎng)格無(wú)關(guān)性的程度后，容差因子對(duì)數(shù)值模擬時(shí)間的影響接近一致。然而研究過(guò)程中也發(fā)現(xiàn)過(guò)高的容差因子不會(huì)縮短數(shù)值模擬時(shí)間，對(duì)于網(wǎng)格數(shù)為14 472的模型，容差因子增加至4 反而增加了數(shù)值模擬時(shí)間，對(duì)于網(wǎng)格數(shù)為18 276的模型，容差因子設(shè)置為4 會(huì)導(dǎo)致計(jì)算過(guò)程不收斂。綜上所述，長(zhǎng)期冷卻計(jì)算過(guò)程中采用的容差因子為2。

圖3 不同網(wǎng)格數(shù)和容差因子數(shù)下的模型中軸溫度分布Fig.3 Temperature distribution on the axis of calculation domain under different mesh number and tolerance factor

3 計(jì)算結(jié)果分析

3.1 速度場(chǎng)和溫度場(chǎng)

圖4和圖5分別為計(jì)算時(shí)長(zhǎng)24 h 和96 h 時(shí)堆芯捕集器內(nèi)的速度場(chǎng)和溫度場(chǎng)，計(jì)算時(shí)氧化物層和金屬層的初始溫度分別設(shè)置為2 400K和2 200K。由于金屬層內(nèi)沒(méi)有衰變熱，其溫度迅速降低900K 以下，并在3 h左右就已經(jīng)完全凝固。氧化物層中上部在衰變熱和冷卻邊界作用下形成了不斷變化的自然對(duì)流，雖然自然對(duì)流的強(qiáng)度不高，但在其攪拌作用下該區(qū)域未發(fā)生明顯的熱分層。由于熔融物黏度隨著溫度降低以指數(shù)形式增加，在氧化物層下部區(qū)域的熔融物溫度較低，因此黏度很大流動(dòng)緩慢，導(dǎo)致熱分層現(xiàn)象十分明顯。

圖4 熔融池內(nèi)速度分布 (a) 24 h，(b) 96 hFig.4 Distribution of velocity field in corium pool (a) 24 h, (b) 96 h

圖5 熔融池內(nèi)溫度分布 (a) 24 h，(b) 96 hFig.5 Distribution of temperature field in corium pool (a) 24 h, (b) 96 h

圖6 為不同計(jì)算時(shí)長(zhǎng)下與氧化物層接觸的堆芯捕集器側(cè)壁面導(dǎo)出的熱流密度。側(cè)壁面受熔融池沖刷的位置熱流密度普遍在0.56MW·m-2以下，隨著側(cè)壁面處熔融物凝固結(jié)殼熱流密度降低至0.2MW·m-2以下，在計(jì)算時(shí)長(zhǎng)2 h 時(shí)，局部最大熱流密度也僅達(dá)到約0.67MW·m-2。坩堝式堆芯捕集器的側(cè)壁面的導(dǎo)出熱流遠(yuǎn)低于大容積沸騰臨界熱流密度，可以認(rèn)為側(cè)壁面不存在燒毀風(fēng)險(xiǎn)，堆芯捕集器的底部壁面與氧化物層之間由金屬層隔離，因此，導(dǎo)熱量很小也不存在燒毀風(fēng)險(xiǎn)。

圖6 不同計(jì)算時(shí)長(zhǎng)下的熔融池平均速度和平均溫度Fig.6Heat flux density of corium pool under different calculation durations

圖7為在計(jì)算時(shí)長(zhǎng)12～96 h的范圍內(nèi)，堆芯捕集器內(nèi)熔融池的平均速度和平均溫度的變化。隨著計(jì)算時(shí)間增加，熔融池平均速度呈現(xiàn)波動(dòng)降低的趨勢(shì)，其中的波動(dòng)過(guò)程可能是氧化物層內(nèi)的自然對(duì)流不斷變化導(dǎo)致的。在當(dāng)前計(jì)算時(shí)間范圍內(nèi)，衰變熱由約13.3MW下降至約7.8MW，而捕集器內(nèi)平均溫度下降速率由12 h 時(shí)約5K·h-1迅速降低至24 h 時(shí)的約2K·h-1，到96 h 時(shí)已經(jīng)降低至1.1K·h-1左右，可見(jiàn)當(dāng)前事故緩解策略的降溫過(guò)程十分緩慢。在事故末期氧化物層完全凝固后，系統(tǒng)的熱量導(dǎo)出能力將更低，因此采用坩堝式堆芯捕集器作為嚴(yán)重事故緩解策略時(shí)，需要相關(guān)設(shè)備和支持輔助系統(tǒng)在很長(zhǎng)時(shí)間內(nèi)保持可運(yùn)行性。

圖7 熔融池平均速度和平均溫度Fig.7 Average velocity and average temperature of corium pool

3.2 熔融池結(jié)殼形式

圖8 為計(jì)算時(shí)長(zhǎng)4 h、12 h、24 h和96 h時(shí)熔融池內(nèi)的固相分布。在氧化物層冷卻過(guò)程中，上部、側(cè)壁和下部的結(jié)殼厚度均逐漸增加，同時(shí)熔融池內(nèi)的最小固相分?jǐn)?shù)也會(huì)逐漸增加。

圖8 熔融池內(nèi)固相體積分?jǐn)?shù) (a) 4 h，(b) 12 h，(c) 24 h，(d) 96 hFig.8 Soild volume fraction of corium in corium pool (a) 4 h, (b) 12 h, (c) 24 h, (d) 96 h

值得注意的是，在氧化物層上部被冷卻的大部分熔融物在完全凝固之前會(huì)受重力影響并隨著自然對(duì)流向底部轉(zhuǎn)移，這導(dǎo)致上部的結(jié)殼厚度的增加不大，而高固相份額的熔融物聚集在氧化物層下部導(dǎo)致該區(qū)域的結(jié)殼厚度顯著增加。側(cè)壁上凝固的熔融物也體現(xiàn)出向下遷移趨勢(shì)，導(dǎo)致側(cè)壁上形成了上薄下厚的結(jié)殼形式。氧化物層內(nèi)硬殼呈現(xiàn)上述形態(tài)是因?yàn)楫?dāng)前熔融物固液相線間的相變溫差很大，遇冷后不會(huì)迅速完全凝固，因此存在受重力和流動(dòng)影響而遷移的現(xiàn)象。在當(dāng)前結(jié)殼方式影響下，氧化物層向側(cè)壁的導(dǎo)出熱量占與總釋熱率的比值，由12 h 時(shí)的約60%逐漸降低至96 h時(shí)的約41%?？紤]到頂部與水直接接觸部分的換熱效果更好，這種結(jié)殼方式對(duì)提升熔融池的冷卻速率是有利的。

4 結(jié)語(yǔ)

采用多物理場(chǎng)耦合計(jì)算軟件COMSOL 搭建了VVER 堆芯捕集器計(jì)算模型，對(duì)捕集器中熔融池的冷卻過(guò)程進(jìn)行了數(shù)值模擬研究，主要結(jié)論如下：

1）在熔融池分層結(jié)構(gòu)下，金屬層會(huì)迅速凝固，氧化物層中上部會(huì)形成不斷變化的自然對(duì)流流動(dòng)且溫度分布較為均勻，氧化物層下部流動(dòng)緩慢并出現(xiàn)了明顯的熱分層。

2）在氧化物層上部被冷卻的大部分熔融物在完全凝固之前會(huì)受重力和自然對(duì)流影響向下部轉(zhuǎn)移，導(dǎo)致上部結(jié)殼厚度增加緩慢，而底部結(jié)殼厚度增加迅速。

3）坩堝式堆芯捕集器的安全裕量較高，冷卻壁面不存在燒毀風(fēng)險(xiǎn)，但其中的熔融物冷卻過(guò)程十分緩慢，需要相關(guān)設(shè)備和支持輔助系統(tǒng)在很長(zhǎng)時(shí)間內(nèi)保持可運(yùn)行性。

作者貢獻(xiàn)聲明朱光昱負(fù)責(zé)數(shù)值模擬和文章撰寫(xiě)；郭超、劉巧鳳負(fù)責(zé)文章撰寫(xiě)；李春、依巖負(fù)責(zé)工作支持。